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AUTOMATION AND COMMUNISM 


Translated from Avtomatika i Telemekhanika, Vol. 22, No. 11, 
pp. I-IV, November, 1961 


On October 31, 1961, the work of the 22nd session of the Communist Party of the Soviet Union concluded — 
a session dealing with the establishment of communism, which opened a new chapter in world history. Communism 
is already close to fulfillment in the near future, and this session opened the door wide not only for our descendants, but also 
for the present generation of the Soviet people. In the new program of the CPSU, there is proclaimed for the first 
time a scientifically-based, concrete plan for the formation of a communist society, which will ensure world-wide 
Peace, Freedom, Equality, Fellowship and Happiness for all nations. 


In his speech on the Program of the CPSU at its 22nd session, N. S. Khrushchev said “In the next two decades 
in the USSR, we will create the material and technical basis for communism, The solution of this main economic 
problem will be the basis of the main line of our party." 


In studying the historical documents of the 22ndsession of the CPSU, we must be guided by the assertion of 
V. 1. Lenin concerning the necessity of transforming communism from prepared, learned formulas, advice, pre- 
scriptions, regulations, and programs, into its main aim — that of uniting our direct work and of converting com- 
munism into a guide for our practical use. 


Soviet scientists must be led in their work by this most important problem posed by the 22nd session; “In 
the development of Soviet science, we must strive for a level of excellence that will make it possible for us to 
attain the leading position in all main directions of world science and technology..." Among these directions, 
of special importance is that of automation, which displays the highest promise of the ability to increase the 
productivity of labor. 


In the conditions under which communism is being created, automation opens up a new era in the development 
of machine technology. Machines for complex mechanization and automation form, in the words of Marx, the bone 
and muscle of the communist production system. 


Automation will bring about a massive growth in the productivity of labor not only because it frees man from 
the necessity of directly controlling each individual machine and process, Thanks to automation, the development 
of production and science (the deciding factor in the powerful growth of the productive forces of society) will not be 
limited by the capabilities of human beings, whose rapidity and accuracy of reaction, memory, and capacity for 
taking into account multiple conditions limit, in many cases, the possibilities of scientific and technological pro- 
gress, Automation opens the way for the development of new processes that can be carried out with high speed and 
great precision, 


In the resolution of the 22nd session of the CPSU, according to the report by the Central Committee, it was 
noted that soviet scientists have achieved great success in the region of cybernetics, and in the development of 
automation and telemechanics, It was also pointed out that a complex mechanization and automation of the pro- 
cesses of production were necessary. 


In the new CPSU program, the aim is proposed of carrying out, in the next twenty years, the large scale, 
complex automation of production, with a more rapid transition to automatic factories and enterprises, in order to 
ensure a high technical and economic development. There must be an acceleration of the adoption of high-per- 
formance systems of automatic control. We must achieve broad application of cybernetics and of computing, 
decision-making, and control devices in the production processes of industry, in the construction industry and in 
transport, in scientific investigations, and in the calculations involved in the accounting and statistics necessary 
for planning and for projecting new construction. It was stressed in the CPSU program, that the development of 
theoretical investigations in the realm of the complex mechanization and automation of production, as one of the 
definite regions of technical progress, leads to a number of the most important problems in science. 


1285 








At the present stage — the stage of creating the material and technological basis for communist society — 
purely empirical methods not solidly founded on theory, can no longer be used to solve the large scale scientific 
problems implicit in the CPSU program, We must now commence the fundamental investigations in automation 
that will ensure that certain ultimate goals will be reached, such as the complex automation of production which 
in the future will become complete, the creation of controlled thermonuclear reactions, the mastery of cosmic 
space, etc. The main direction of theoretical development must be shifted towards the branches of investigation 
that deal with fundamentals. Contemporary control theory must not be limited to the research into and the de- 
velopment of those systems already in existence. One of its main purposes must be the inculcation in practical 
workers of a broad outlook, leading them to methods of seeking new ways of constructing technical automation 
systems and of creating methods by synthesis. 


Control theory must, to a large extent, be used in profound scientific exploration and for the clarification 
of what machines can do, what they cannot do, and how we can widen their range of application. In the task of 
constructing highly effective automatic systems, the development of technical cybernetics is of great importance, 
since the latter is the theoretical basis for a large complex of technical disciplines — automation, telemechanics, 
electronics, computing techniques, etc. Technical cybernetics is the most highly developed part of cybernetics; 
its methods contribute to the promotion of science, not only in the technological realm, but also in economics and 
biology. 


At the present time, when, as the resolution of the 22nd session of the CPSU has stated, the Soviet Union is 
entering the epoch of development of scientifically-based, automatically-controlled production processes, the 
central problem is that of optimal control. It is therefore necessary to develop the theory of optimal automatic- 
control systems, and to relate this theory to the problem of the construction of automatic-control devices. 


Many results already exist concerning the theory of optimal systems, these results having been obtained by 
applying the calculus of variations. Such general principles and methods have been derived in this way, for example, 
the maximum principle and dynamic programming methods, The application of these methods to concrete pro- 
blems has, however, not yet emerged from the initial stage of development, and some problems have been found 
to be so complex that only a general approach to their solution has been so far attempted. 


Another important direction is that of the development of the theory and the principle of operation of systems 
containing automatic devices, as applied in cases where the conditions are variable. This class of systems contains, 
for example, automatic filters with variable characteristics, which are self-correcting for variable conditions, and 
which ensure the best possible separation of a signal from a background of random interference. 


The principles and theory of devices for the automatic investigation of objects and the recognition of form 
must be developed. Such devices, to some extent, will be modeled on the processes occurring in man's nervous 
system during the recognition of form. 


The theory must be developed of self- teaching, self-adjusting, self-organizing systems, in which only the 
first steps have so far been taken. 


To the important theoretical directions to be investigated, we must add that of the development of complex 
automatic-control systems, such as multi-coupled systems, in which the control must be organized so that, in the 
presence of complex intersecting couplings, several controlled parameters satisfy prescribed conditions. 


The development of automatic systems for the control of logical operations requires a further advance in the 
theory of finite automata — a theory in which a clarification will be made of the possibility,in principle,of con- 
structing devices of discrete action, and in which methods will be developed by the analysis and synthesis of 
systems that operate according to prescribed operational algorithms. Such a theory will make it possible to build 
automatic systems that will function in a way that approximates the functioning of human beings in the sphere of 
intellectual activity. 


Human beings will take part in many complex control systems. The functions of the automatic controls and 
of the human taking part in the control process will be so closely intertwined that it will become impossible to pose 
and solve problems of control by the consideration of the technical aspects of these problems only, without taking 
into account psycho!ogical and physiological factors that are relevant in the man's part of the control process, It 
will be necessary, to ensure the maximum effectiveness of the joint operation of the man and the machine, to take 
into consideration the limitations of both the man and the machine. Here the man must be employed under con- 
ditions such that his work will be productive and interesting, but not exhausting. Much work will have to be carried 
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out on the creation of machines that will analyze the characteristics of the objects of production, machines that 

are self-adjusting, and machine-automata, in order to lighten the work of humans in calculations, All these devices 
must be developed toward a range of greater complexity of application and a higher level of automatic operation in 
planning. 


Our era will see a rapid growth of knowledge. A literal chain reaction is evident in the growth of knowledge 
at the present. Serious problems are arising in connection with the automation of the storage and transmission of 
vast quantities of information. The work of humans in looking up information and in its processing must also 
be greatly simplified, and this work will be partly carried out by specialized machines called automatons. Such 
machines will be a powerful aid to scientific progress, and will form a new branch of the national economy — that 
of the scientific research industry. The development of automatic-control theory involves several branches of 
mathematics, such as variational methods (classical and direct), methods for the approximate solution of problems, 
statistical methods from probability theory and mathematical statistics, the theory of games, etc. Theoretic - 
logical methods must also be used — logical algebras, predictive calculus, and also various methods from functional 
analysis. For solving complex practical problems, developments are needed in approximate methods of calculation 
and in methods of programming the solutions of problems for mathematical computing machines. This all indicates 


the necessity for a close liaison with mathematicians, both pure mathematicians and specialists in numerical analysis 
and programming. 


Automation will not be possible without progress in the technical components. In this region, the key problem 
is that of improving reliability, and this problem must be solved both by developing more reliable components and 
methods of combining these components, and by seeking methods of constructing reliable systems from unreliable 
components. A wide application will be necessary of semiconductors, ferrites and ferroelectric materials, thin 
crystalline films, gas-filled tubes, etc., for developing miniature, high-reliability automation devices, A great 
perspective opens up here for the development of quantum physics. Thus, the main optical range of radio waves 
can be used to construct controlling devices with extremely large flows of information. The solution of problems 
of this type will demand the organization of close cooperation between physicists and specialists in radio electronics. 


The problems of carrying out the new program of the CPSU require a broadening of the scientific research 
front, both with reference to the problems of automation of the more distant future (two decades distant), and with 
reference to problems having solutions that will produce great economic results in the near future. To attain these 
results, it will be necessary to broaden and improve the network of scientific and experimental- constructional organiza- 
tions developing new principles for the construction of automatic-control systems and new cybernetic devices and 
automation systems, which is under the control of the central establishments for the guidance of public construction 
and the Councils of National Economy. This also applies to the network of scientific laboratories and institutes con- 
nected with heavy-industry establishments. Scientific research into automation must also be given a more im- 
portant place in the institutions of higher education. Decisive measures must be taken for the strengthening and 
perfecting of the material base of our knowledge of automatic control, and for the attraction of the most capable 
and creative forces into scientific activity in this field. 


The instrument-manufacturing and machine-manufacturing industries must be oriented toward the creation 
of single state system of automation and telemechanical devices, which is necessary for the contemporary complex 
problems of automation. This system must be developed on the mass-production principle. A relatively small 
number of different items will yield the possibility of organizing mass production, will decrease the cost of the 
individual items involved in automation, and will raise their reliability. This will permit the acceleration of the 
technical retooling of the national economy. 


It is of great importance to ensure that conditions are such that there will be a rapid introduction into practice 
of new scientific results. In view of the new problems, we cannot be satisfied with the existing rate of introduction 
of the present methods of automation and telemechanics into industry. Up to the present, there has been no standard 
form of organization of scientific research and of projected work on automation that would guarantee any real 
coordination between the work of the scientific, the experimental and constructional, and the planning organizations, 


We must organize as rapidly and as widely as possible the practical coordination of the creative scientists, 
engineering technologists, and planners for the development of new techniques, in particular in connection with the 
complex problem of drawing up an outline of future projects to be undertaken by the communist society in basic 
industries. This must be done with the aim of building a more completely automated industrial complex, and of 
introducing automation in the sphere of organization of control and cooperation of industrial activity. 
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Soviet scientists and engineers working in the area of automation are proud of the part they have played in 
the establishment of a high rate of automated production under the communist regime, and in the creation of con- 
ditions for the complete development of the human personality and the harmonious relation between separate in- 
dividuals in our whole society. 


Soviet specialists in automation warmly support the splendid plan for communist growth presented by the 22nd 
session of the CPSU and will spare no energy in endeavoring to bring about the realization of the great communist 
ideals. 
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ANALYTICAL CONTROL DESIGN IN SYSTEMS WITH RANDOM 
PROPERTIES 


Ill, OPTIMUM CONTROL IN LINEAR SYSTEMS. MINIMUM MEAN-SQUARE ERROR 


N. N. Krasovskii and E. A. Lidskii 


(Sverdlovsk) 

Translated from Avtomatika i Telemekhanika, Vol. 22, No. 11, 
pp. 1425-1431, November, 1961 

Original article submitted March 18, 1961 


We consider the problem of determining the optimum control law under the assumption 
that the control system is described by linear equations. The minimization of the mathematical 
expectation of the mean-square error is taken as the quality criterion for transient processes. The 
solution of the problem is based on the formation of the optimal Lyapunov function, as a quadratic 
form with coefficients that are functions of the random variable. The question of the solvability 
of the problem is discussed. 


1. The Problem* 


1. In the present article, we will assume that an automatic control system is described by the linear equations 
of the disturbed motion 


dz 
= Mu (M) +... + ain(M) en t+ (ME + T(z, &) (i 
E=§(%,-.+, In» n)- 


ere (1.1) 


Here x; (t) are the coordinates of the system, and the variable n(t) describes the random properties of the 
system and is a Markov process, 


The interference y,; is assumed to be a pulse random function with frequency ). The dependence y; = 
vis (Vi (t), py (x, &) of the random magnitude of the interference on the mismatch x(t) and the controlling 
signal €(t) (the control law) is assumed to be known. The quantities yj are random with known dispersions 
6? = M {v2} and the coefficients kj; define the statistical relation between yj and yj (M (vyvj} = 1436;6)). 
For the linear system (1.1), we take 


i = Pit +--+ + pinTn + Pind (4; = const). (1.2) 


2. The quality of a transient process is measured by the size of integral 


I (Xp, n= Sae{] > ©P + EP (a4, . » os Ins 1) / %y = Tyo, +++ In = Ino, N= No 


i=1 


for ¢ = ft, = o}} dt. (1.3) 


The problem is to find a control law € = €°(x, n) that will ensure the probability asymptotic stability of the 
motion x = 0 for the Eqs. (1.1), and the minimization of the integral (1.3). 


* We will use the results and notation of [1, 2]. 
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2. Method of Solution 


1. We will seek the optimal Lyapunov function ({1], section 4; [2], section 1) in the form 


= 2) bis (n) 2423. 


i, jt 

In order to satisfy the conditions that must hold for the functions v’, we must require that the coefficients 
bij) satisfy, uniformly with respect to n, the Sylvester criteria [3] for the sign-definiteness of quadratic forms, 
and that the following relations hold: 











(% )., + 3 x3 fe E02 _ 0, (2.1) 
(4 ae me 3! ae + poe __ min in | os » a + | : (2,2) 
i=] =i 


2. We now write Eqs. (2.1) and (2.2), assuming that n (t) varies within finite limits 1, << y(t) <A», and 
using the formulas given in (2) (Eqs. (1.5) - (1.7): 





am) [eves + Youn a5] AE” +\ v¢e, Badan 0 


i=1 j=1 1 


Ov (x, n) 
—q(n)vo(z, +z >y “Fade, Pibts5i65- (2.3) 
he A 
We consider two particular cases. 
A. The quantity n(t) takes a finite number of values 4y,... nm, with transition probability IPs | (1), 


section 2), In this case the argument 7 will not be written, and the functions will be provided with the editaigt l 
when n =n 
1° 


B. The function q (ny, 6) has the density p(n, 6). 
The interference y;, and y j Will be assumed to be independent for i = j, in which case the Eqs. (2,3) become: 





for the case A 
4 ve “ ° 
t=) 2 aa me py ats | 
j=1 
1 ~~ e 
mel ix1 i 
for the case B 





n Ay 
“T= Des E m+ Yay (m 2 + »(@, B)p(n, B)dB 


i=1 j=1 


ho 
—9(n)v(*, n) +> 2 = pict. (2.5) 


The equations for v’ and ¢° are obtained from (2.4) and (2.5), taking (1.2), (2.1), and (2.2) into account. 
For example, for the case B: 


1290 








me; 


2.5) 








Ae 
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it j=1 ; 
n n 2 n 
. 2. 
+525 | Hints’ + > wits | 0? + >) a+b =0, (2.6) 
fon j=1 i=1 
aF oy av dty = aw a 


3. The Calculation of the Coefficients of the Form v® (x, n) 


We consider case A. We substitute the quadratic form into the Eqs. (2.1) and (2.2) and eliminate €. We thus 
obtain the relation 


n n n k n 
23) |( 3 ols) (D alls) ]+ 3 rm SS Of — a) a2 
i=1 mel 


>) wi) | 


j=1 , 


n n n n 2 
| > (4 > Hie, ) +A > Cone > rei) | 


i=1 j=1 i=1 j=1 


+>) E Dat ( 
i=1 





= + > H=0. 
i+ p> Oe aa 1 


i=1 


When we compere coefficients of identical products, we obtain a system of 1/7, kn(n + 1) algebraic equations 
(i a ’ 
for the coefficients bij (2=1,..., kos 


an k 
D OB al? + of af) + Ab ripest) + Dd) Pay, (Ss? — bY) - 





i=—1 ml 
n n 
[Sot nae Hints) | Pa Minti 
i= {a= 
! = (3.1) 
i+’ > OS a 
i= 
“+ om (p, = 1,..., W). 
0 (p+#s) 


For the case B, we obtain a system of integral equations for the functions bin) (these equations can be 
obtained from (3.1) by replacing the summation with respect to m by an integral with respect to 6). The equations 
obtained can be generalized for the stochastic case to yield the equations obtained by A, M. Letov (4), part IV). 

4. An Approximate Method of Solution 


1. The method is described in (2) ection 2). We will illustrate this method as applied to the linear system 
(1.1). We consider the auxiliary system of equations 
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oi = [iy () Ty 4 eee -f din (n) Ly] 6 + BEC; (n) +V 8%; (x, §1) 
+ (1—#)§& —(1 — 8) 2; 
(i=1,..., n) (E=&: for @ =‘). (4.1) 


The parameter 9 is introduced in such a way, that for 9 = 1 the system (4.1) becomes the system (1.1), and 
for 9 = 0 it degenerates into the simple equations * 


Ft — my + & for #=0. (4.2) 


We will try to find the control ;,...,€, that minimizes the integral 


1 (8) = (mM D2? + B+ (1 — O) DET / 21 = Boy. «Zn = TaN = No 


for t = tft (4.3) 


n 


For (4.2) the optimal function v’ is easily obtained in the form >) ,,2?. For 9 > 0, we will seek the 
i=1 





function v° as a form with coefficients depending on 9 (except for the dependence on n= ny,..-.7 , » The 
derivative dM {v} /dt, in view of (4.2),has the form 
a Pe si 32 Ov nt eps +(1—#8)(§; —2z) +9 Days 
inn” j=1 
+4 5 Pim (Um — ¥1) + +9 S 2 re 202, (4.4) 
me¢l i=1 Ox 


2. We will solve the problem by a known method (2), section 1) (with the special feature that the con- 
trolling signals €,,...,€, must be eliminated after the partial differentiation with respect to all the variables 
&4 of the equations corresponding, in the present case, to (2.1) (see the footnote). We obtain the equation 


2 > (3 2 bs bia : [© cos —(i— 9%) x |} 
[orn 


k n 
+05) Pm dD) OP — dl?) xa; + 00 > 


ns i, j=1 i=1 — 
n n n n 2 
' l > 
[e 5! (« > Hs, +(1—) DPirar+ Or Y, (sisting 2 wi) 
i=1 jim i=1 =1 =1 





1+0A >) (OPS? Hin 41) 


i=1 


—(1 —@) >) ( > Hes) + >) a3 = 0. (4.5) 
i=1 


i=2 ‘ j=1 


* The process of control for $ = 0 breaks up into n independent channels, each of which has its regulator. When 

® varies from 0 to 1, a continuous change in the system takes place, and for 9 = 1 control is obtained by a single 

action €, for all coordinates x;. The solution in the case of n controlling actions is similar to that described for a 
min 


single action, with the difference that the sum corresponding to (2.2)has € for all €. 
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A further differentiation of (4.5) with respect to 9 is carried out, the coefficients of identical products xix; 
are equated, and the resulting equations are solved for the derivatives @b{{)/00. The result is the system of 
dif :rential equations 

avs" ql) (m) (m) : 5 
36 = > j rr ore (i==1,..., a; f == 1; m=, ...h). (4.6) 

These eqns must be integrated with respect to 9 over the interval 0 <9 <1, taking into account the 
initial values bt (9 =0). The {pjtial conditions (wy 0)=0 for ¢j, bi) — V2—1) are obtained from Eqs. 
(4.5) for 9 =0. The solutions b;, (9) for 9 = 1 yield the coefficients of the optimal Lyapunov function for the 
original problem. The question of the possibility of integrating Eqs. (4.6) is related to that of the existence of a 
solution of the optimal problem. This question is considered in the following paragraph. 


5. The Existence of an Optimal Solution 


1. The question first arises of the conditions under which in probability it is possible to ensure the stability 
of the prescribed motion x (t) = 0, by the choice of the control law € (an attainable control). After the solution of 
this problem, we can consider the problem of the existence of an optimal ¢°, minimizing the mean-square error. 
We will not discuss the question of the existence of an attainable control here. A further special article will be 
devoted to this question. 


2. We will show here that, with the assumption of the existence of an attainable control, we can prove that 
there exists an optimal control. We note, first of all, that the existence of an attainable control ensures the con- 
vergence of the integral of the mean square error, if the system is linear (5). The basis of the possibility of con- 
structing an optimal control, when it is known that an attainable control exists, rests on the above-described method 
of constructing auxiliary systems with parameter 9. We give the corresponding reasoning, assuming that the inter- 
ference depends only on the mismatch in (4.1). The optimal function v (x,, 9) for the system (4.1) must satisfy 
the equation 
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Here the expressions for £, and €, are 


n 
1 | Ov . Ov 
eS fe ae peas | ed ' (1) OU; 
n= — 4[Ha—o+0 Bag], 
: 1 dv, . 
FF . +o for ae 


We differentiate (5.1) with respect to 9, set 8v,;/ 89 =a, and transfer to the right-hand side the terms not 
containing a, or the derivatives da; / 0x; : 
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Using (4.4) and the formulas for €, and €;, we note that the left side of this equation is the derivative of the 
mathematical expectation dM {a,} /dt, obtained from (4,1), and the right side is some quadratic form, Thus, (5.2) 
can be written in the form 





os e)) = s 8) (®) aj (L=1,..., ). (5.3) 


i, j=1 


On the basis of the results in [5], we may now assert that, if the system (4.1) is asymptotically stable in the 
mea.’ for the values we are considering, then for an arbitrary quadratic form on the right side of (5.3) there exists 
a unig te solution — the quadratic form a; - of the equation (5.3), i.e., we can calculate from Eqs. (5.3) the 
derivat. ves and hence the derivatives 


Th possibility of continuing the solution v (x, n , 9), known for 9 =0, to all values of 9 > 0 up to 9 =1, 
is now proven in the following way. The reasoning above shows that the derivatives 0b; /09 exist and can be 
calculated from Eqs. (5.3) for all values of 9 for which the optimal problem is solvable and there exists an optimal 
Lyapunov function v’, The values of these derivatives 0b;; / 0% can then be used to continue the solution to larger 
values of $, This method of continuation of the solution v (x,n , 9) can be carried out only up to the first value 
9° at which either it becomes impossible to solve Eq. (5.3), or some of the coefficients bjj (9) ~ © for9 > 9°. 
It becomes impossible to solve (5.3) for 0b;; / 09 only at the approach to a value 9° such that the optimal system 
loses its asymptotic stability in the mean, i.e., when the function v for 9 +9* ceases to be positive definite. This 
however, is impossible, since v for 3 < 9° is ‘bounded below by a a positive number at the points 2x}= 1 because 
of the fact that v is equal to the eget I (x9, 1 9, 9), which possesses this property. It is also impossible for 
bj; ~ @ for § + 9°, since for 9 < 9° we have v=1, (x, 7 9, 9), and this quantity is bounded above at every point 
Xo, because of the existence of a permissible solution for which lL, < Iner- We therefore see that the process of con- 
tinuing the solution of v (x, n , 9) up to 9 =1 by the integration of the equations for 0;; /0% is possible. This 
proves the existence of an optimal Lyapunov function, and hence an optimal control law for the original problem. 
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The paper studies the theory of formulating difference equations for stationary and non- 
stationary discrete shaping filters. A specific example of the application of discrete shaping filters 
in the statistical analysis of the nonsteady-state movement of a pulse system is studied. 


1, Statement of the Problem Involving the Determination of a Discrete Shaping 


Filter 


In many problems involving the statistical investigation of pulse systems it becomes necessary to obtain 
random discrete processes with specified correlation functions. 


A random discrete process with a specified correlation function can be obtained by means of a linear pulse 
system. In fact, assume that a pulse system has the operator Ay. Then the output variable Y of the system when 
its input is subjected to a perturbation X is the result of applying the operator A, to the function X: 


Y = A,X. (1) 


If X is a random function (in the general case it is a nonstationary random function) with the correlation function 
K,([n; m], then the correlation function K,[n; m] of the random function Y is a result of the double application of 
the operator A to the correlation function K, (once with respect to the variable n, and the second time with respect 
to the variable m): 


K, [n; m] = ApAmKg[n; m]. (2) 


Formula (2) demonstrates that if the correlation function K, [n; m) is specified, then the correlation function 
n; m] can be obtained by an appropriate selection of the system operator A, Thus the random process Y with a 
specified correlation function can be obtained from any random function X. 


However, the simplest form of the operator A arises in the case where the random process Y is derived from 
discrete “white” noise that consists of a discrete random function with uncorrelated values.° 


A pulse system that converts discrete “white noise into a discrete random process with a specified correlation 
function shall be called a discrete shaping filter. 


Assuming that in the general case a discrete shaping filter can be called nonstationary, we shall write its 
equation in the form 


®,, (A, n) Y (n] = ¥n(A, n)V [n}, (3) 


where V is discrete “white” noise with a unit intensity; Y is the “shaped® random process with the correlation function 
Ky (n; m); © and ¥ are difference operators which are defined by the relationships 


* We study a class of random processes such that they can be obtained from a solution of difference equations. 
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k h 
®,, (An) = - a{nj Ai, WV,(A, n) = NV bj [nn] A). (4) 


i=0 i-0 


where at is the operator which governs taking the difference of order i with respect to the argument n; aj and bj 
are real functions of the discrete argument. 


Thus, if the correlation function K, [n; m] of the shaped random process Y is specified, it follows that the 
problem of determining the shaping filter reduces to determining the difference operators @p(A, 0) and wp(A, n). 
We shall proceed to solve this problem. 

2. The Solution of the Problem 


The difference equation for the stationary shaping filter (or its frequency response) can be found using formulas 
for the expansion of the spectral density of a discrete random process which were derived in [1]. However, in the non- 
stationary case this method cannot be applied. Here, we shall cite a general method for determining a discrete 
shaping filter; the method is applicable both in the stationary and nonstationary cases. 


Since V is discrete “white noise" from which the process Y with the desired correlation function Ky (n; m] is 
formed, it is possible to say (see Appendix) that Ky [n; m] satisfies the following difference equations for n > m: 

@, (A, n) K,,[n; m] = 0 (n>m), (5) 

@,, (A, m) K,[n; m}] = V,, (A, m) VW, [n; m] (n>), (6) 


where w/,(A, m) is a difference operator which is conjugate with Wp (A, n) and is determined by the relationship 


h 


Win (A, m) W,[n; m) = > ( -1)'ai, {b,[(m — ij W, [n; m— i}}, (7) 


i=0 


and W, (n; m) is a function which is determined by the solution of the homogeneous equation 
@, (A, n)Wy[n; mj} -=0 (n>m) (8) 
for the initial conditions 


AW, in; m| =0 (=O, i,..., k—2), AN W,, [n; m| es CRT (9) 


n=m n=m ay [m] . 


Equations (5) and (6) are the original equations for determining the operators @ and ¥. 
Since Ky [n; m] is a solution of Eq. (5) for the corresponding initial conditions, it is obvious that 
k 


K,(n; m| == »y ni lm) Ay ind, (10) 


{==} 


where ), is the fundamental system of solutions for (5), nj are functions which are determined by the initial con- 
ditions, k is the order of the difference operator ©. 


We shall first determine the operator @, (A, 0). It is possible to prove that for linear difference equations the 
following theorem is valid. 


Theorem. The fundamental system of solutions );[n) (i=1, 2,.. . , k) fully defines the linear homogeneous 
difference equations of k-th order with a coefficient equal to unity for the highest-order difference. This equation 
is written as follows: 
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A, [n} Ag(n}... Ax in] y [n) 
AA, [n] Adgin]... Ady(n] Ay{n] 


2 46..¢ © 2 IRs 6 £2 2, C7 ee oe eee 


A* a, {n} — A* Ag inj... A™ Ay in} A* yfn} 


W, (A, n)y[n} = = 0. (11) 


From formula (10) and the theorem formulated above it follows that the difference operator @ , (A, n) can be 
determined if the specified correlation function K, [n; m) is represented in the form (10), and the determinant of 
the type (11) is then formulated and expanded with respect to the elements of the last column. It is obvious that the 
order of the difference operator  p(A, n) is equal to the number of functions ) ; into which the specified correlation 
function for the shaped random process is expanded. 


We shall now proceed to determine the difference operator ¥,(A, n). Based on (8), (9) and (11), it is possible 
to demonstrate that the function W, [n; m) is defined in terms of the fundamental system of solutions ), in the 
following manner: 








A, [m] Ae [m]... Ae [m] 
(—1)k—1 AA, [m] Ads [m]... Ady [m} 
W,[n; m] = acta all ed sn? boalit ts fe eee ’ (12) 
A*-*),[{m} A*—*A,{m]... A**Ay [mn] 
A, [n] Ao [mn]... Ax [n} 
where a), is the coefficient of the highest-order difference in the operator 4, 
A, [m + 1) Agim-+ 1)... An (m + 1] 
AA, [m + 1] Ad,g{m-+ i]... Ad, [m + 4] 
DM] 4 oss Lies « segues be Xe agetibets be (13) 
A*—7,{(m +4] A*7A,{m+1]... A Ag tm +4] 


Knowing the difference operator @, (A, n), it is possible to find the result obtained by applying the product 
of the conjugate operators W,,(A, m) and ¥ 4, (A, m) to the function W, (n; mj in accordance with formula (6). 
Then,determining the form of the product of the conjugate operators ¥ and ¥° and expanding it into conjugate 
multipliers, we obtain the sought-for operator Wp, (A, n). 


The difference operators ® ,(A, 0) and ¥ yp (A, n) which are found in this fashion exhaust the solution of the 
problem involving the determination of a discrete shaping filter. 


As an example, we shall find the difference equation for a filter which forms the oft-encountered nonstationary 
discrete random process with a correlation function of the form 


K, [n; m] = @ [rn] @ Im] jee. (14) 


from “white* noise which is of unit intensity. 
For n > m we obviously obtain 
a [n] = @ [n] e~*", m1 [m] = @ [m] e®™. (15) 
Formulating a determinant of the form (11) and expanding it with respect to elements in the second column, 
we obtain (we omit the common multiplier e~ “") the operator which determines the left side of the difference 
equation for the filter: 
®,, (A, n) y [n] = @ [n] Ay [n] — @ [n + 1] €~* —@ [n) y [a]. (16) 


In order to determine ¥, (A, n), we have 
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a, [m] = @ {m], D|[m] =@ |m+ 1] eam +), 


Therefore, 


Wo ln; m] = @ (n] —a(n—m—1) 
° pim]@(m +i] ° 





and in accordance with formula (6), we obtain 
¥, (A, n) v [n] = @ [n] @ [m] (Av [n] + v [m)). (17) 
On the basis of (16) and (17) it is possible to write the difference equation for the sought-for filter in the form 


@ InJy([n + 1] —@ [n + 1] y [n] &* =@ [n] @ [m] v [n + 4]. (18) 


Equation (18) is essentially an algorithm in accordance with which we obtain a random process with a correla- 
tion function of the form (14) from the uncorrelated sequence V. 


We shall now study a case involving the application of discrete shaping filters. 


3. The Use of Discrete Shaping Filters for the Statistical Analysis of Nonsteady- 
State Motion of a Stationary Pulse System 


The problem of the statistical analysis of pulse systems reduces, as a rule, to determining the dispersion (or 
more rarely, the correlation function) of the output variable of the system when its input is subjected to a random 
process. Here we shal! present a method for determining the dispersion of the output variable of the stationary system 
in a nonsteady state mode. 


Assume that a stationary pulse system with a weighting function W([n, ¢€] is specified in which the transient 
responses end after a time t = NT, where T, is the period required for closure of the pulse element switch and N is 
an integral positive number. For the system under study it is required to determine the dispersion Dy {[n, €) 
(0< n+ e< NT,) of the output variable Y of the system when its input is subjected to the random process X 
with a correlation function K,[n; m). 


The weighting function of a stationary system is a function solely of the instant t=n+ € at which the signal 
is sampled at the output. However, in the time interval 0 < t< NT, (i.e , during the transient response) the system 
should be treated as a nonstationary system. In general, a stationary system that is at rest until the instant t = 0 and 
is excited at the instant t = 0, can be treated as a particular case of a nonstationary system whose parameters vary 
stepwise at the instant t = 0, 


Studying lattice functions [2] for the sake of simplicity, we can demonstrate [3] that the dispersion of the 
variable of the nonstationary system is determined from the formula 


D, (nj = > W[n; m,] > W [n; m,] Kx [m,: mo). (19) 


™m,=—0o ™m,y=—0o 
In the particular case where X is discrete “white” noise with the correlation function 
K,,.[n— m] = D,s([n — m], (20) 


where the o-function is determined by the relationship 


for n=, 


1 
sin—mi =| 4 for n=, 


the expression for the dispersion becomes 
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D,|(nj=D. SN W?{n; my). (21) 


Thus, if the input signal of the system is a random signal of the “white® noise type,then the dispersion of the 
output variable is determined very simply. 


Making use of discrete shaping filters, it is possible to replace the investigated system whose input perturbation 
is a random process having the correlation functions K,[n; m] with an equivalent system whose input perturbation is 
discrete “white* noise, Here the equivalent system is a serial connection of a discrete shaping filter with the 
weighting function W,[n; m] and the original investigated system. 


In that case the dispersion of the output variable will be computed from formula (21), where the weighting 
function W(n; mJ of the original system should be replaced by the weighting function We,[n; m) of the equivalent 
system. 


Formula (21) is valid for determining the dispersion for either a stationary or a nonstationary system. Here we 
shall cite a different method which is applicable only to systems with constant parameters in the case where the 
input perturbation is stationary and has a correlation function of the form* 


Kx [n— m] = e-7!"—™ qa. (22) 


Applying the method cited above, it is possible to show that a filter which shapes the random process X with 
the correlation function (22) from “white noise” V is described by the equation 


AX [n] -|- (1 —e~*) X [n] == V1 — eV [n + 1). (23) 


The weighting function for such a filter is determined by the expression 


Wy in—m| = Vie ee) 


(n> m). (24) 
Thus,the original investigated system with the weighting function 

W can be replaced with an equivalent system formed by a serial con- 

vi Pn dO Twine | 'y nection of a filter with the weighting function (24) and the original 


“lm 1 alae system (Fig. 1). 


Weq [n,¢] 


Under these conditions, we must keep in mind the fact that the 
random process X is obtained from discrete “white” noise V in accord- 
ance with the formula 

bf [rn] = N' Wig in — mV [m}. 


m-=—oo 


Therefore it is possible to study the case where the inputs of the system may arrive before the instant n = 0 
(i.e., before the instant at which the system is turned on; m < 0), or after the system has been turned on (m > 0). 
This fact is indicated in Fig. 1 by the presence of a switch which closes at n = 0. 


For a signal o [n - m) at the filter input,the output variable of the filter will be determined from the relation- 
ship 
s[n—m| = V 1 —ee- z(in—™) (in > m). (25) 


The weighting function Wed: m] for the over-all system will evidently be determined as follows: 
Wein: m) = Y: Win —m,]z[m, — m) =e Vi-—e 2 >» Win—m,je72"" (m< 0), (98) 
m=) men 


* An analogous method for determining the dispersion for continuous systems was developed by G. H. Laning (4). 
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since the signal at the input of the investigated signal with the weighting function W begins to arrive from the in- 
stant at which the switch closes (n = 0). 


For m > 0 the input signal of the investigated system begins to arrive from the instant at which the o-function 
is applied rather than from the instant at which the switch closes. Therefore, we shall have 


n 
Veg [n; m) — Vi —— @-os »s Ww [n ly, m,) e—2 (Mm—m) (m > 0). (27) 


m=m 


From expression (26) it follows that the function 


F [n| = > W [n — m,]| e~?™ (28) 


m,=0 


is the response of the system to the signal e~" for n > 0 (the system is turned on), and to a zero signal for n < 0. 
Thus expression (26) can be written in the following form when (28) is taken into account: 


Weqin; m}] = Vi —e~%e2M™F [nn] (m <0). (29) 


For m > 0 we obtain 


n—m 


Weqin; m) = Vi — ee > W [n — m — m,] e~2™ = V1 —e-** F[n—m!| 


m,=0 


(m> 0) (30) 
from (27). 
Therefore, we have 
—22 Lam 
Weqin; m] = (Vs vied een F [n] (— oc <n g—}), 
) i_¢e™ F [n — m} (O< magn). (31) 


Substituting expression (31) into formula (21) for the dispersion, we find 


oj] 


n 
Dyin} = >) (4—e-*) 2” F2 [n} + D>) (1 —e-22) F?2 [Jn — mJ. (32) 


m=—oo m=0 


Performing the summing operation in (32), we obtain 


Dy, (n| = e-2*F? [nj + (1 — e-2) >> F? [m]. (33) 


m=0 


Since F[n] is the response of the investigated system to a signal of the forme" “” for n> 0, it follows that 
the dispersion D,[n] can be defined as a time function using the model whose block diagram is shown in Fig. 2. 
The method for determining the dispersion (for € = 0) reduces to the following under these conditions. 


A signal of the form e~%" is applied to the input of the 
model for the pulse system when n=0. The response of the 
model which can be represented by a function F[n, ¢] is applied 
to the clamping unit @ whose interrupter operates synchronously 
and in phase with the switch in the pulse element of the system. 





Fig. 2. 


The output quantity of the clamper is a step function 
whose values over the interval n <t <n+1 are equal to the values F[n, €] at discrete instants of time (« = 0). 
The step function F[n, 0] is applied to the square-law device K and then to the storage device £, the scale blocks 
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and the adder. It is evident that the output quantity of the entire model (i.e,, the dispersion) is a function of the 
discrete argument n=0,1,...). For n = N a steady-state movement of the system occurs and the dispersion will | 
be constant in this case. 


The dispersion of the output variable at the instant t = n+ € can be obtained from the formula 


Dy (n, &| = e-*F? [n, e| + (1 —e-*) 5) F?[m, el. (34) 


m= 


In order to determine D,[n, €] for € = €, (0 <€ 9 <1) it is possible to use the model described above,with the sole 
difference that the interrupter of the clamper must not operate in phase with the switch in the pulse element of the 
system but must close the circuit with a time lag € 9. 


The model shown in Fig. 2 can be represented in a different form if we note that the input e"™" to the model 
is the weighting function of a filter described by the equation 


Ay {n| 4+- (1 —e-*) y[n] = z[n + 1]. (35) 


Thus, by substituting the response of this filter to a o- function for the input e~ “" and changing the sequence 
in which the filter and the system are switched on (this is possible since their parameters are constant), we obtain 
the model shown in Fig. 3. 





Fig. 3. 


The block diagrams which have been cited for the models which determine the dispersion of the output variable 
of a stationary pulse system in a transient mode are general ones. They are applicable both in the case where digital 
computers are used and in the case where analog computers are used. In the case where digital computers are used 
the corresponding blocks in the diagram are networks for the numerical solution of the equations that describe the 
operation of the blocks. 


The method cited above for determining the dispersion can easily be extended to the case of arbitrary 
stationary perturbation, since it is easy to show that the correlation function K,[n - m] for a stationary process can 
be represented in the form 


! 
K, [n = m| - > A. at | n—m | (36) 


¥=1 
with an arbitrarily high accuracy. 


In conclusion the author expresses his sincere appreciation to Ya. S. Itskhoki, A. P. Grishin and V, N. Gotesman 
for the attention which they devoted to this paper and the help which they gave him in its preparation. 


APPENDIX 
Assume that a nonstationary pulse system is described by a single difference equation of the form 
®,, (A, 2) y-[n,e] == ¥,, (A, 0 & In, 0) (0<e < 1), (37) 
in such a way that the weighting function W(n, «; m] for this system satisfies the difference equation 


Y, (A, n) W [n,e; mj) = ¥, (At) o [n — m, 0] (0g e < 1), (38) 
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where the step function o is determined from the relationship 


1 (mgnqgm-+ t), 
o [n — m, 0] =| (39) 
0 (m>n>m-+ 1) 


It can be demonstrated that the correlation function of the random process Y at the output of such a system, 
when its input is subjected to a nonstationary random perturbation X with a correlation function K, [n; m],is deter- 
mined from the formula 


n ™m 
K, In,e; me] = >) Win,e,t) >) W Im,es; a1 K, (i, 0; a, 0). (40) 


l=—oo A=—0co 


Assume now that X is discrete *white" noise with an intensity identically equal to unity. In that case ex- 
pression (40) can be written as 





n 
| >) Wiae; 1 WImei tl) (n<m), 
ie testes | ona (41) 
>) W (ne; W (mer; 1] (n> m). 
l=—oo 


Applying the operator ,, (A, n) to the second Eq. (41) in the variable n, we obtain 


n 
®,(A, n) Kyin,e; m,e:) = YS) ®,"(A, n) W [n,e; 1) W Imes; 1] (n> m), 


l=—oo 


whence in view of (38) and (39), we have 
m mm 


Dd W ime; 1) ©, (A, n) W (ne; = >) W Im,er, NY, (A, ) o [n—m, 0]=0. 


{=—oo l=—oo 


Therefore for n > m the correlation function for the random process at the output of the pulse system satisfies 
the following homogeneous difference equation when the input signal is discrete “white” noise: 


®,, (A, n) K, [n,e; m,e:) = 0 (n> m) (42) 


for the corresponding initial conditions, 


From the first of Eq. (41) it follows that Ky {n, €; m, €3] for n < m represents the response of the system to a 
perturbation that coincides with the weighting function of the corresponding conjugate system. Thus, we have 


®,, (A, n) K, [n,e; m,e1] = ¥,, (A, 2) W [m,ei; 2] (n<m). (43) 


Treating K,(n, €; m, €,) as a function of m, it is possible to demonstrate in analogous fashion that it satisfies 
the following nce equations: 


;, (A, m) Ky (n,€; m,e1] = 0 (m>n), (44) 
©, (A, m) K, (n,¢; m,e1) = ¥), (A, 0) W [n,e;m)  (m< a), (45) 


where T=m+ €, and 6° and y® are difference operators that are conjugate with the operators @ and ¥, respectively. 
For € = €, = 0 formulas (42) were used in solving the problem of determining a discrete shaping filter. 
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METHODS FOR REALIZING OPTIMAL FILTERS WITH 
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(Moscow) 
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pp. 1441-1452, November, 1961 

Original article submitted April 24, 1961 


The paper studies methods for realizing an optimal filter with a finite memory which was 
defined in the Zadeh-Ragazzini problem. It is demonstrated that discrete filters can be used to 
achieve a better approximation of an optimal system than the approximation which can be 
achieved using filters with an infinite memory. The proposed method for realizing filters with a 
finite memory is illustrated using the example of the realization of an optimal filter- differentiator. 


1. Statement of the Problem 


In (1) - [2] a study is made of the problems involved in the synthesis of an optimal (in the sense of the mean- 
square criterion) system which is designed for isolating the useful signal g(t) from noise n(t) and transforming g (t) 
in accordance with a specified operator L(p). It is assumed that n (t) is a stationary random time function with a 
fractionally rational spectral density, and that the useful signal is an unknown regular time function which can be 
expanded into a Taylor series 





i 
g(t —t) = Does ty (1) 


In (1) it was assumed that the number of terms in the series (1) is limited. In that case, the optimal system 
assures a minimum noise dispersion at the system output (oSpr) and a zero average value of the error involved in 
reproducing the useful signal L (p) g (t): 


Ope min, M (8 out opt“)! = L(p) g(t), 


where Bout opt is the useful signal at the output of the optimal system. 


The pulse transient function h,,,, (t) for the optimal system applies to the class of filters with a finite memory 
and is a solution of the following integral equation when (r +1) conditions are satisfied: 


T 
| R(t — gg(t) dt = 0+ Tt + + + Te”, (2 
0 
T 
Ho = \ Pop dt) dt = Xo, 
T 
Hy we thoy f*) dt = — %, 
ns 
= "Th, = xX, (— 4 (3) 
p, | wag (ree %, (— 1) 
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where R(t— r) is the correlation function for the noise; y 9, y4,..., y, are indefinite Lagrange multipliers; yy 
is the i-th moment of h,.(r)and «4 =0 or xj =1 depending on the form of the operator L(p) (for example, for 
Up) = const we have %)= 1; x= x, = ...=%,= 0); T is the magnitude of the memory of the optimal filter and 
ris the power of the highest degree term in the series (1). 


The optimal pulse transient function is equal to [1]: 


Noplt) = Ao + Art +--+. + A,t' + C8 (t) + C28! (t) +... 4 Cyd (0) 
+ D,b(t — T) + DW (t —T) +... + Day, (t — T) 
+ Bye +...4+ Baye?! for O0<t< T, 
hop ft) =0 for t<0, t>T, (4) 


where Ag... A,; By... Boy; Cy...Cir4y; Dy... Dy4, ate coefficients whose values are determined for a 
specified R(t - r) by substituting expressions (4) into Eqs. (2) and (3). 


In [2] a general case was studied in which the number of terms in the series (1) was unlimited. It is demon- 
strated that when an optimal system with h,,,, (t) determined according to expression (4) is used, the average value 
of the error with which the signal is reproduced is no longer equal to zero: 


M [Zour opr (t)] =L (p) 8 (t) a Spr (!), 
where € op, (t) is the dynamic error with which the useful signal is reproduced. 


The same paper estimates the magnitude of the maximum dynamic error €opt max- In [2] - [3] it is demon- 
strated that for an appropriate choice of the magnitude T of the filter memory it is possible to assure the required 
relationship between Sopt and € opt max (i.e., it is possible to assure a correspondence between the noise stability 
of the filter with a finite memory and the degree to which the useful signal is distorted. In [3] a method is given 
for approximating the series (1) with a series that has a finite number of terms (a method for choosing the value of 


eZ. 


The practical realization of an optimal filter with a finite memory is difficult. From expression (4) it is 
evident that for the simplest form of hop: (t) (C, = Ca =... = Cyy, = Di = ... = Dyy,= 0) the corre- 
sponding optimal transient function for a filter with finite memory consists of a combination of pure lag sections 
and integrating sections that are combined by means of multiplier and adder units (3). We know that integrating 
sections, and especially pure lag sections, are in themselves complex and difficult to realize. Moreover, a practical 
system will be stable only when the coefficients of the individual terms and the expressions for the optimal transfer 


function correspond exactly to the coefficients of the optimal transfer function and all the mathematical operations 
are performed without error. 


In this paper we study methods for using the characteristics of a filter with a finite memory to approximate 
the dynamic characteristics of analog and digital engineering elements that have a comparatively simple practical 
realization. Based on the proposed method we demonstrate the possibility of realizing an optimal filter-differentiator 
(L(p) = T,, pl. 
P 


2. Approximating the Dynamic Characteristics of a Filter With an Infinite Memory 
Using the Dynamic Characteristics of a Filter With a Finite Memory 


In the general case we can assume that stable linear systems with a pulse transient function h, (t) of the 
following form have a simple practical realization: 


ha(t) = 1 (t) >) ((boi + Suit +... + but!) cos wat 


i=1 


+ (cot + cuit +... + cnt)! sin ot] ec! + GB (2), 
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for t< 0 


a iso M7 (5) 


0 
1(t) =}, 


It is necessary to establish the conditions which must be used to make hy (t) approximate ho», (t). 


Usually a function of the form h,,,,,(t) is approximated by a function of the form hp (t) on the basis of the 


opt 
mean-square criterion [4] 


\t hye) —h,(t)/? dt = min. (6) 
0 


It is known [4] that by increasing the number of terms in the series m and the degree of the polynomials | 
in expression (5) it is possible to achieve any degree of accuracy for the approximation of ho,,;(t). However, in 
practice we cannot complicate the approximating filter without limit (m and J are limited). Therefore, the useful 
signal L(p)g(t) at the output of such a filter may be appreciably distorted. For the form of g(t) under study, the 
dynamic error €,,(t) at the output of the approximating filter depends on the moments of the pulse transient function 
Ry (t). 


First, we shall study the case where g(t) is represented in the form of a Taylor series with a finite number of 
terms. In order to limit the value of €,(t) the first r moments of the pulse transient function for the approximating 
filter must satisfy conditions analogous to conditions s (3), However, the moment hy(t) (a filter with infinite memory) 
will be functions of time over the entire interval 0 <t <a. Therefore, the additional conditions of the form (3) 
are satisfied for only one value t = ty: 


te 


Pao li, = \ h,(t) dt = %, 


to 


Mni | —t, = \ tih,(t) dt = x(—1)', 


ty 


Por los, = \ Th, (t) dt = x,(— 1)’, (7) 


where pj is the i-th moment of h,, (r). 


Taking into account the fact that the limiting value of yj corresponds to t, = o, it is expedient to assume 
t,=a. In contrast to filters with a finite memory,only the limiting value of €p(t) (for t-— a) will equal zero in 
the case under study, and for all other values of t > 0 the value of €p(t) will be nonzero and will depend on the 
first r moments. 


If the series (1) is not limited, then ¢€,(t) also depends on moments of a higher order; here no conditions at all 
are imposed on the moments beginning with the moment of (r+ 1) order. For both forms of g(t) the dynamic error 
is determined by the expression 


e q(t) = L(p) g(t) —\ 2 (t— 1) h(x) ae. (8) 


In order to assure correspondence between the dynamic error and the noise stability of the approximating 
filter (just as for the optimal filter (3)) it is expedient to introduce the condition 


| ®n max | = Sp (9) 
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where €p max is the maximum value of the dynamic error at the output of the approximating filter, and op is the 
mean-square value of the noise at the output of the approximating filter. 


The introduction of conditions (7) and (9) makes it possible to isolate those systems which will reproduce 
L(p) g(t) with the minimum dynamic error from the class of filters (5) and for which the noise stability will corre- 
spond to the degree of distortion that appears in the useful signal. If hopr(t) is approximated on the basis of the 
mean-square criterion using the approximating filters indicated above, 4 for bounded values of m and | the 
approximating filter will reproduce L(p)g (t) with the minimum dynamic error which can arise at the output of 
stable filters with an infinite memory. 


Conditions (7) and (9) can be treated as a system of equations relative to the parameters h,(t), Therefore the 
approximating filter must have a minimum of (r + 2) arbitrary parameters. 


3. Approximating the Dynamic Characteristics of a Discrete Filter by Means of the 
Dynamic Characteristics of a Filter with a Finite Memory 


Assume that a certain signal y (t) = g(t) + n(t) is applied to the input of an optimal filter with a finite memory. 
The signal at the output of the filter with the finite memory is yop¢(t) and is equal to 


t 


y (t)= \ hk. (ti—t)y(t)dr. (10) 
opt \ opt 


Assume that we are interested in the values of Yopt‘t) at discrete instants of time under conditions where the 
difference between the values of the argument (the time t ) corresponding to two adjacent values of Yop is equal 
to the magnitude of the filter memory T: 

itd 
Yope(ti+1) = \ Peal tits —t) y(t) dr. (11) 


K 
On the right side of expression (11) we obtain a definite integral instead of an indefinite integral; this appreciably 
simplifies the realization of such systems. 
We rewrite (11) in the form 


T 
Yop tits) = \ Moe DY (ti4, —t) dr. (12) 


0 


In the integrand of (12) the mirror image of y (t) relative to '/, (tj4, — ti) appears instead of the function y (t); 
i.e., the absolute magnitude of the argument (the time t ) decreases, and this is inconvenient for purposes of practical 
realization. It is evident that the magnitude of the integral appearing on the right side of expression (12) is not 
altered if we take the mirror image of the weighting function h,,»(r) relative to YT instead of the mirror image 
of the function y(t) and shift the origin of y(t) by the amount T in the direction of a decreasing argument (we take 
ti4, instead of t ): 


T 
Yo C+) = a —t)y(t& +t) dr. (13) 


Expression (13) will be made the basis for simulating a discrete filter with a finite memory. Such a system, 
irrespective of the form of Hope 7), will consist of the following elements: 


1) An integrator which performs the operation of integrating a certain time function over a time 0 <t <T; 
2) A device which integrates the mirror image of the weighting function over a time 0 =t <T; 
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3) A multiplying device in which the continuous multiplicacion of two functions is performed over a time 
0st xT. 


It follows from expression (4) that the weighting function h,,,,(t) can be represented as the sum of a certain 
continuous function (Bop e(0) over the interval (0, T) and a §-function with its derivatives. 


Obtaining h,,,,(t) [the mirror image of Dope 0] reduces in the general case to reproducing a certain time 
function according to a specified law over the interval [0, T]. 


Of the remaining components of the weighting functions only the 6 -function is comparatively easy to realize 
(more precisely, the integral of the 6 - function): 


T 
AYont= \ [C15 (T— t) + Did (x) y(t + t) dt = Cry (ts + 7) +Dry (ti). 
i) 


Thus, it is possible to realize any discrete filter with a finite memory with an accuracy of up to derivatives 
of the §-function. Individual sections of such a filter (an integrator, a device for reproducing AY opt» etc.) can in 
principle be realized either on an analog basis or on the basis of digital techniques. 


We should underline the fact that in contrast to an optimal filter with a finite memory we are dealing with a 
digital system (the output signal changes at discreze instants of time). The errors in the operation of the individual 
elements of such a digital system can affect only the accuracy with which a specified operation is performed; they 
cannot affect the operational stability of the system (here errors are not integrated). At the same time the dynamic 
error in reproducing the useful signal L(p)g(t) may increase in comparison to the error of an optimal analog system. 
The dynamic error ¢ q(t) at the output of the discrete filter, in contrast to € opti) [2], is in the general case equal 
to the maximum value of the difference 


T 
&q (t) = L(p) g (t) —{ Rls (tits — t) dt... 
for tip <t<ty,+ 7. (14) 


For example, for L(p) = const, the value of € opt) is equal to 


Pry 
e —— | o(r+1) 
le IIS aT lets | 


in accordance with (2); here g“*) is the maximum value of the (r+ 1) -th derivative of g(t). 


In determining the maximum value of ¢ q(t) we assume that during the time interval [tj+,, t{4,+T]) the 
(x + 1)-th derivative remains equal to the maximum value (the worst case): 


Pets | (r-+1) 
le (Ol< (7 + “) ewan, 

The additional dynamic error at the output of a discrete filter with a finite memory can be reduced in two 
ways. 


1, By designing the filter in the form of a complex unit that includes several discrete filters of the type 
described above. These discrete filters operate with a “phase shift" relative to each other. The additional dynamic 
error is reduced by a factor k, where k is the number of discrete filters. 


2. When a single discrete filter is used, the magnitude of its memory T must be chosen while taking into 
account the additional dynamic error [expression (4)]. As computation shows, under these conditions the additional 
dynamic error can be reduced by approximately 50%. 
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In practice, it is necessary to estimate the expedient accuracy with which hop, (T - 1) is reproduced and 
Ay opt {8 measured when a discrete filter with a finite memory is realized. For this purpose, it is necessary to com- 
pare the magnitude of the dynamic errors and noise stabilities of a discrete filter with a finite memory for ideal 
and approximate reproduction of h(T - t) (and for measurement of Ag,,,,). 


4. Realization of the Optimal Filter-Differentiator 


The operator L (p) of the optimal filter- differentiator is equal to L(p) = T,»p. 


We shall study the approximation of an optimal filter- differentiator corresponding to a noise signal with the 
correlation function by the dynamic characteristics of analog and digital differentiators. In that 
case the optimal pulse transient function of the filter-differentiator is equal to [3] 


6 ’ 
hope (t) = Fara ear pin |a?T — 20% + (2+ aT) d (t) 
—(2+4+aT)d(t—T)]) O<t<7), 
Rope (t) =0 (t<0, t>T7). (15) 





A. Approximating h op; q(t) by means of the pulse transient function hy, ,,(t) of a filter-differentiator with an 
infinite memory. The dynamic error at the output of the filter-differentiator with an infinite memory is determined 
practically by the zero and second moments of the pulse transient function for the system* in the case of a slowly 
varying useful signal: 


(0) 
lenl S18 Bool +18 Hnal for T<t < om, (16) 


where gt) is the maximum value of the useful signal, gi?) is the maximum value of the second derivative of 


the useful signal, pp» ts the zero moment of h,, po and pn? is the second moment of hp piv. 


Table 1 cites the expressions for the pulse transient function, the first moment, and the dynamic error for the 
filter- differentiator with a finite memory in which the numerator of the transfer function is equal to T,p and the 
denominator is equal to a polynomial in p of the first and second degree. 


The filter- differentiator of the first order with an infinite memory has two unknown parameters: the co- 
efficient of the derivative T,, and one root of the characteristic equation. The second-order filter-differentiator 
correspondingly has three unknown parameters: T,, and two roots of the characteristic equation. 

The coefficient T,, of the pulse transient function cited in Table 1 is determined from condition (7). 

As computations have shown, of the three forms of h, ,(t) for the second-order system,the best approximation 
of Dont q(t) can be obtained using the pulse transient function of a second-order system for multiple roots of the 
characteristic equation. The unknown parameter ay (a root of the characteristic equation) for first and second-order 
systems can be derived from condition (9). 


Assume that Eq. (9) is written in terms of the maximum value of the dynamic error over the interval [T, oo]. 
In accordance with conventional methods of analysis, it is necessary to find tp max for which | &,|/ “ye”, reaches 
a maximum (for the condition that T < tn max < ©). From the expressions |e, |/|g‘?),,|, cited in Table 1 for 
first and second-order systems (with multiple roots) it is evident that in both cases tp max will be determined by a 
transcendental equation. Therefore & max | / |\g‘?),,| can be found only for specified numerical values of 
Sinax!/| &(2),,| and oy using approximate methods (for example, graphically). 

The following sequence of solving Eq. (9) should be noted, We substitute the value of o,/8 from Table 2 
(which cites the noise dispersion at the output for two types of filter-differentiators with infinite memory) and the 
corresponding value of |@n|/| 2 (from Table 1) into Eq. (9). For the assumed numerical values of the co- 
efficients a, \emay /\g2,,| and B/|e{).,.| of Eq. (9) we used the method of successive approximations to de- 
termine the value of oy satisfying Eq. (9). In finding the specified value of a4,the values of op are found using the 
formulas in Table 2, and the values of | €,, ,,,| ate determined graphically (as indicated above). 


* Expression (16) is valid for 0 <t <, But in order to compare |¢,,| with | €ope| we shall determine |¢€,,| over 
the interval T <t <@. 
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If we take into account the fact that the functions &n max (4, |e. \/| se... \)and % (a, B/| Sinax !) 
remain continuous within bounded limits when the arguments a, |g), | / /\g@.. | and B/|g@),.| vary (for any 
a), it follows that from the particular values of a, determined above rit is possible to find the general laws govern- 
a the variation of «a, for a continuous variation of these arguments. Figures 1a and 2a show the functions a, = 

(a) for constant values of B/|g®, | and |g, |/|g@..| (|g. |/|e@.,|=9 and | earl /18 Tax 

=e. Oo Nes max max 
B/|g.x|=1 and B/|gun = 10) 



































































































































TABLE 2 
T f Relative magnitude of the noi 
yee Gee dispersion at the output (0,/8°) 
; 2 
first order ren 
filter-differentiator a Oy 
of the analog type : 
second order Lets 
(multiple roots) Verto) 
—at 4 
first type 2 fee 
discrete filter- d! 
differentiators stg 
second type i late —(i-e @)I 
d? 
of oop 
i a 
=] oe am Ig? '— Fa 
i tee 
12 | Se tS! 6b ae. 
4 | 7 
Wao tt] Spams 
|g? 1 ig - 1 ORE as a8 
Wyner] eae) ST eta 
6 max StH st 92s 
V-Faat Co AAA | | 
4 woo — <a 9 max 4 
2h ae 6 ol Waa at 
T's ot ta far ee ee 
a b 
o? Jt 
Fig. 1. The functions — = @n(o) and ay Fig. 2. The functions ioe = 9p (a) and ay 
¥ nla) for — re oT] I. = ¥n(a) for -——— 5 Te] 
of tee SS ee ee *— represents the first order; 


—-------- represents the second order; 
represents the optimal system. 


ee ee oe represents the second order; 
represents the optimal system, 








For the same values of B/|8,,! and |g@,.1/|8,,| Figs. 1b and 2b show the curves for the relative 
value of notes dispersion 6? /B? = ?, “(a). at the output of the e filter-differentiator with an infinite memory. These 
values of oF, 2/8” were computed from the formulas cited in Table 2; the values of ay were determined from the 
curves in Figs. 1a and 2a. In Figs. 1b and 2b, we have also shown the curves which characterize the noise stability 
of the optimal filter- differentiator. 
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From these graphs it is evident that the noise stability of a first-order filter- differentiator deteriorates in com- 
parison to the noise stability of an optimal filter-differentiator when a is increased or when | cael \g coal is 
increased, For the same values of od [\g? and B/| a the noise stability of the second-order filter- 


5 
differentiator proves to increase in comparison to the first-order system as a increases, For a — oo the noise dis- 
persion at the output of the first-order system increases without limit, while the noise dispersion at the output of the 
second-order system remains a bounded quantity. From Figs. 1b and 2b it is evident that the noise stability of the 
second-order filter- differentiator is close to the noise stability of an optimal filter-differentiator for | 8°) . | / 
|giax! —+0. The case where this ratio is much greater than 1 [this is characteristic for a slowly-varying 
function g (t)]is of great practical interest. 


B. The approximation of B opr,  (t) by the pulse transient functions h,, q(t) of discrete differentiators. Table 1 
cites the expressions for the pulse transient function, the first moment, and the dynamic error for two types of discrete 
differentiators, In [5] it is demonstrated how the practical realizations of these systems can be achieved and their 
dynamic characteristics are determined. 


Discrete filter-differentiators of the first and second types have two unknown parameters (Tp, 7 gy and Tp 
T qa» respectively), The coefficients of the derivative for both types of discrete differentiators is chosen on the basis 
of condition (7). The first moments of the discrete differentiators (Table 1) do not depend on time for t > T. The 
coefficients T,, of the pulse transient functions for the discrete differentiators cited in Table 1 are determined from 
condition (7) which is satisfied for any t > T. 


The unknown parameters rq, (a discrete differentiator of the first type) and rq» (a discrete differentiator of 
the second type) are determined from condition (9). 


a 
14 


12 
l 
08 









\ 











04 
G...F 4 $§$:+6¢ a 
b a 
Fig. 3. The functions oy = gq (a) and Fig. 4. The functions oe = gq (a) and 
Tqa= ¥q (a) for (2) Ss Td= ¥q (a) for —5—- =10. -_*—-¢*- 
Sax | | ®max| 
represents a discrete differentiator of the first represents a discrete differentiator of the first 
type; — — — represents a discrete differentia- type; --—-— represents a discrete differentiator 
tor of the second type; represents the of the second type; represents the optimal 
optimal system, system. 


The dynamic error at the output of the discrete differentiator for t > T is independent of time and depends * 
basically on the second moment (the zero moment of the discrete differentiator is equal to zero fort > T). Table 1 
cites the expressions for the dynamic error, and Table 2 cites the expressions for the relative magnitude of the noise 
dispersion at the output of discrete differentiators of the first and second types. 


Substituting the values of | €gmax| and og into Eq, (9), we obtain transcendental equations in rg, (or rq»). 
The parameters rq, and rqp are computed in a manner that is analogous to the computation of oy. Figures 3a and 
4a show the functions ty = tpy(a), and Figs. 3b and 4b show the functions 53/8" = 9, (4) for two types of 
discrete differentiators in the cases B/| gi. .| = 1 and B/|g® | = 10. 


The same figures show the corresponding graphs for the optimal filter- differentiator. 
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It is easy to show that for a =<7.5 the noise stability of an optimal filter- differentiator approximately corre- 
sponds to the noise stability of a discrete filter- differentiator of the first type; for a = 7.5 it corresponds to the noise 
stability of a discrete differentiator of the second type, This value of the degree of noise correlation (a = 7,5) can 
be assumed critical in choosing the type of discrete differentiator for purposes of approximating h,,, ,(t) by means 
of its pulse transient function, When §/|g'T),,| increases, the critical value of a decreases (for B / |g‘), | = 10 
it is equal to approximately 2.5; see Fig. 4b). As computation shows, when a varies within the limits 0< a <10 
and both types of discrete differentiators are used,the noise stability of the approximating system will differ from 
the noise stability of the optimal system by no more than 10 to 15%, 


C, Comparison of the noise stabilities of analog and discrete filter-differentiators. From a comparison of the 
gtaphs in Figs, 1 and 2 with the graphs in Figs. 3 and 4, it is evident that in practice the noise stability of discrete 


filter- differentiators is no less than three times greater than the noise stability of filter-differentiators with an in- 
finite memory. 


For the purposes of illustration, Fig. 5 shows comparative oscillograms for the operation of a discrete filter- 
differentiator of the second type and an analog filter-differentiator (a second-order system), An identical signal 
was simultaneously applied to the inputs of these filter-differentiators; the signal consisted of a useful signal plus 
noise [g (t)+n(t)]. The over-all input signal was oscillographed simultaneously with the output signals from the 
discrete filter-differentiator (Agg2+ Ang») and from the analog filter-differentiator [g,.,,(0) + n,,;(0)]. 


PUIOV OH. ok, 


ee ; 
- Jud *,, Ad Grmax 90085 a 
Vy, Imax ).25v sec 


ya 
pore ee fee 


Of + An 
- a 











b B-0.55¥ =005 eb 
g ft+n (0) Imax ,00085 —— 
out out yt, =5.25 v 89° 


a = \ ks 1 beeen een he ae ioe 
A +a, Git)+ntt) 


8d2 





Fig. 5. Oscillograms which characterize the noise stability of analog and discrete filter- 
differentiators. a) Optimal setting of both filter-differentiators; b) for the same setting 
parameters but for double the mean-square value of the noise at the inputs. 


Oscillogram a (Fig. 5) corresponds to the optimal setting of both differentiators for 8 = 0.36 v and g* 
= 0.00085 v/sec®, Oscillogram b corresponds to the same settings when the mean-square value of the noise at the 
inputs of the filter-differentiators is increased to 8 = 0.55 v. From oscillogram a it is evident that for the optimal 
setting the output signal from the discrete filter- differentiator reproduces the variation in the slope of the useful 
signal with a sufficient accuracy, 


SUMMARY 


1, Inapproximating the weighting function of a filter with a finite memory with the weighting function of a 
filter with an infinite memory on the basis of the mean-square criterion, it is expedient to introduce certain 
additional conditions which must be satisfied by the pulse transient function of the approximating system, On the 
basis of these conditions it is possible to choose from the class of filters with an infinite memory those systems 
which will reproduce a specified useful signal L(p) g(t) with a minimum dynamic error and will have a noise 
stability that corresponds to the degree to which the useful signal is distorted, 


2. The weighting function of a filter with a finite memory can be approximated with an accuracy of up to 
derivatives of the §-function using digital systems, Such a discrete filter will be stable under all conditions (this 
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is in contrast to an analog filter), and the additional dynamic error in reproducing L(p)g(t) which arises due to the 
discrete nature of the variation in the output signal can be reduced substantially, 


3. The noise stability of a discrete filter-differentiator whose parameters were determined on the basis of the 
proposed method for typical characteristics of the useful signal and the noise, is no less than three times as great as 
the noise stability of a filter-differentiator with an infinite memory (for optimal settings), 
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THE STATISTICAL INVESTIGATION OF EXTRAPOLATION 
EXTREMAL“-CONTROL SYSTEMS FOR AN OBJECT WITH 
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We investigate the operation of an extrapolation system of extremal contro] which is designed 
to maintain an extreme value of the signal at the output of an inertialess object with a para- 
bolic characteristic. An approximate expression is found for the steady mean error in tracking 
the extremum, in the presence of random interference at the object output and a drift in its 
characteristic. From a comparison of results obtained under identical conditions, using the 
extrapolation and the gradient methods, a region is found where the extrapolation and gradient 
methods of search for an extremum can be rationally applied. 


INTRODUCTION 





In the majority of the systems of extremal control known at the present, the search for the extremum is 
carried out without using any information that may exist concerning the characteristics of the object of control, 
At the same time, in many cases, the form of the equations describing the given object of control, i.e., the equations 
of the hypersurface on which the desired extremum is sought, can be stipulated beforehand on the basis of known 
theoretical or empirical relations. The unknowns in this case are the constants, or the rather slowly-varying quan- 
tities, occurring as coefficients in the equation for the object. We should note that, even if the form of the equation 
of the object is not exactly known beforehand, in the majority of cases it can be represented with sufficient accuracy 
by a finite power series with unknown coefficients. As soon as the form of the equation for the object is known, 
then we can determine the values of the coefficients occurring in this equation from some experimentally obtained 
data, and then determine analytically the action of the control that will yield an extremum. This method can also 
be used to obtain the essentials of the extrapolation method of extremal control [1]. Methods of finding experi- 
mentally the unknown coefficients in the object equation are given in several articles, for example [2, 3]. 


The question naturally arises of the extent to which supplementary information concerning the object can 
be used in improving the search process and in maintaining the extremum, i.e , the question of the expediency of 
applying extrapolation extremal controls. 


Exact values of the coefficients in the equation for the object, and the consequent exact extremum location, 
can be found only in the case when the object parameters do not change with time and when the experiments with 
the object are not subject to random interference which does not satisfy the given equation. Only in this case do 
extrapolation systems have obvious advantages over step extremal systems. It is evident that the quality of operation 
of an extrapolation system will be worsened as the level of random interference increases and as the rate of para- 
meter drift increases. Under these conditions, the more “delicate® step extremal-control systems can, in many 
cases, guarantee a more accurate tracking of the extremum. On the basis of this doubly theoretical comparison of 
step and extrapolation systems, the opinion has been established (see [4], for example) that it is not practicable to 


apply the extrapolational search method in the presence of any significant level of interference in the object of 
control. 


In the present work, we investigate the influence of random interference and object- parameter drift on the 
simplest extrapolational system of extremal control in the case when the object of control is characterized by a 
quadratic. The results are compared with the characteristics of a gradient system, investigated for similar conditions 
by A, A. Fel’dbaum in [5], and it is found that there is a region where both the above types of system are applicable 
in the control of the class of objects considered. 
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Extrapolation Method of Search with a Constant Test Step 





We consider an extremal control system (Fig. 1), consisting of an object O and a controlling device Y. The 
output coordinate w of the object is a random function of the time, depending on the input variable x: 


w(t) = y(t, 2) + dw (2), (1) 
where § w(t) is a random, stationary, centered function of the time; 
y(t, x) = A(t) [x +1 ())}? + ¢ (t) (2) 


and A (t), 2 (t), and c(t) are time-dependent parameters. 


We also assume that the character of the disturbances due to the drift of 
Tou the parameters is such that the relation (2) can be written in the form 
ty é 





wwii wy 
tT 


twee awa dé 








y (t, xz) = A(x +1, + v,t)? + cy + vel, (2a) 











within the limits of a sufficiently-large time interval, where A, 1 9, Co, vy and 
y vg are initially unknown constants [constant in the interval where (2a) is valid). 
—— y There are therefore certain unknown, finite limits within which the values of 
2/9 


A, Vy and v2 must lie. Moreover, it is known that A > 0. 























Thus (2a) determines the equations of the parabolas, convex below, 


Fig. 1. which describe the uniformly reciprocating motion in the (x, y) plane. 


The problem of control is to determine and maintain a signal x , such that the quantity y(t, x) is minimized, 
i.e., such that the mathematical expectation of the output coordinate w is minimized. 


Although we assume that the desired change in x can be produced instantaneously, we nevertheless assume 
that the measurement of the size of w requires a certain definite time At. In this connection, we assume that the 
control process takes place at discrete times i, that are integral values of the relative time t=t/ At. Thus, from 
(1) and (2a), we have 


w(i) = A [x (i) +1, + ai]? +c, + Bi + dw (i), (3) 


where a = vi At, B = v2 At. 


The quantity At is assumed to be large enough, compared with the time of damping of the autocorrelation 
function of the random interference 6 w (t), so that there is no correlation between the values of § w (i) for different 
i, 


We note that an object, defined as above, is similar (with the exception of inessential features) to the object 
considered as an example in [5], in which the gradient method ‘of extremal control was investigated. 


We write the relation (2) in the form 
y (i, 2) = Ag (i) 2* + A; (i) x + Ay (i). (4) 


Here the current coordinate of the minimum of (4) is given by 





D ite gaeig 


It is evident that an ideal process for tracking the extremum would be obtained if x varied according to (5). 


The desired current viaues of A,(i) and ALi AL i), ALi] can be obtained approximately by means of the method of 
continuous regression analysis (6). These approximate values are obtained from a system of equations [6] obtained 
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from a certain set of n pairs of values of 2(i —t), w (i —t)(t =1, 2, ..., n<i), obtained from experiments 
with the object at the times i- r: ° 


A,(i)n+ A,(i) >) ci—1) + A, (i) >) 2G — 1) = D} wi—v, 


t=1 t=] t=] 
Ao (i) 5) ti — 1) + AL) > 2G —) 4+ AL) dD 2AGi—D 
t= l1 t=] t=] 


I 


> x(i —t)w(i —T), 


Ap (i) 5) 2 (i—1) + A.) D Pi—9) + A, D ti —) 


t™=I1 t=1 t=1 


= >) 2(i—1)w(i—?). (6) 


t=1 
To simplify the solution of the system (6), the above experiments on the object are carried out using a method 
similar to that described in [2, 3], namely by introducing the supplementary, periodic, search motion 


a for i = 3k, 
0 for i=3k+1, (7) 


x,(i) = 
p —a for i=3k-+ 2, 


where k is the period number of the search (k =0,1,2...), anda is the fixed size of the test step, 
Thus, the signal x (i) arriving at the input of the object will be 


(i) = 2u(i) + 25 (0) (8) 


where x,, (i) is the signal arriving from the control device. 


The change necessary in the control signal x, (i) will be produced at the time i = 3Nj, where N is some fixed 
number, j=0,1, 2,... Therefore at the times i[3NVj <i< 3N (j-+-1)], which we call the j-thsearch cycle, 
the control signal will remain constant and equal to xy (j). The measurements x(i) in the j-th cycle experiments 
will therefore be relative to this quantity x,,(j), which is constant during the cycle, i.e., we have introduced a new 
coordinate X (i), the initial reading of which varies from cycle to cycle. 


Thus for 3Nj <i<3N(j+ 1) we have 
X (i) = z(t) — tu), 


y(t) = A[X@+LW* +e), (9) 
where L(i) = /(i) + zu (/). 


From (3), (7), and (9), for experiments in the j-th search cycle, we obtain 


ws (j,k) = A [a + L(j) + 3ka]* + ¢ (7) + 3B + bu (i, k), 
ws (j,k) = A[L(j) + (3k + 1) a]* + c(/) + (Bk + 1)B + du (j,k), 
ws (j,k) = A[—a + L(j) + (3k + 2) a}* + ¢(j) + (3k + 2) B+ dw4(j, k), (10) 


* The determinant of the system (6) is nonzero, if at least three values of x (i - r) are different from one another, 
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where k=0,1,..., N—1, 
w, (j,k) = w(3Nj + 3k), we (j,k) = w(3Nj + 3k + 1), 
ws (j,k) = w(3Nj + 3k + 2), 
and L(j) and c(j) are the values of L (i) and c (i) for i = 3Nj. 
The value of the ratio in (5) for time i = 3N(j +1) will be written as 
. Ai (i +1) 
1)=-——_ .. (5a) 
z(j + 1) 2, 04-1) 
The quantity z (j + 1) obtained from (5a) will be used as the increment of the control signal for the (j + 1)-th 
search cycle, i.e., 


Su (j + 1) = tu (7) + 2(7 + 1). (11) 


Here the initial reading of the coordinate X (i) is displaced to the point x, (j + 1), and the experiments of the 
(j + 1)-th search cycle will be described by the system (10), where 


L(j +1) =L(j) + 3Na + 2(7 + 1). (12) 
In a similar way, at the time i = 3N(j + 2) the quantity z(j + 2) will be found, the new control signal x, (j + 2) 
= x, +1) +2(j + 2 will be obtained, and a new displacement of the initial reading X (i) will be carried out. 
We obtain the value z (j + 1) from the results of the j-th search-cycle experiments. 
If we now substitute in (6) the value of X(i) from (7) - (9), and use the fact that in the case we are considering 


3N(j+1)—1 3N(j+1)—1 3N(j+1)—1 
y w= YF X=, > X*(/) = 2Na? etc, 
i=3Nj i=gNj i=gNj 


then, we obtain the solution of (6) in the form 
N—1 








p> (wi (7, &) —ws (i, &)) 
aed a 
N—1 
ss >} (wr i, &) + ws (7, k) — 2we (7, kD] 
A; (j + 1) pre —= 2Na? . (14) 


From (14) and (10), after some transformation, we obtain 





AG +1) =¢ (44-2) [24°] 
N-1 
—2B+ 4S} (dw (j, k)—dws(j, XI}, (13a) 
k=0 
N-1 
Hy (j +b) = AS + ahs Fld (7, b) + bus (j, be) — 2ur (i, h)]. (14a) 
k=0 


It is obvious that the tracking of the extremum at the instant i = 3N(j + 1) will be performed more accurately 
when the value of z(j+ 1) approximates more closely the quantity L(j)+ 3Na, i.e., when the values of A, (j + 1) 
and A,(j +1) are obtained more accurately. It can easily be seen from (13a) that there is some finite number N 
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that for given L(j), a, a, 8, and a given level of interference, for which we will obtain the most reliable value of 
A,(j+1). At the same time as this follows from (14a), the most reliable value of A(j +1) is obtained for N - o. 
Noting also that A (j + 1) does not depend on L(j), i.e., it is independent of the strength of the control signal and 

the position of the initial reading of X (i), we deduce that, to obtain the most reliable position of the extremum of 
A, (j +1), it is desirable in the general case to determine from the data the number N, of whole search periods, where 
N, > N; here the choice of N, does not depend on the length of search cycle N determining the period of variation 
of the control signal and the number of experiments for the determination of A, (j +1), In this case 


N-—1 


As (7 +1) = AHF 4 i tbeay 7, ) + bang (7, hk) — 2davy (i, bp. (4b) 


=—(N,—N) 





For sufficiently large N,, the random interference will have practically no effect on the size of A, i+). 


It must however be noted that, for fixed N, (N, > N), the quantity A, (j+1) can be obtained only for j+1 
= — » and this latter condition necessitates either the omission of the control on the original search cycles, 
i.e., the slowing of the search process, or the introduction of a variable value of N, into the first search cycles, i.e., 
a complication of the extremak control device. 


In view of the above discussion, we will limit our investigation to the case N, = N, which is the simplest from 
the point of view of the realization of the theory. 


The search algorithm, in this case, will be obtained from the equations (7), (8), (5a), (11), (13), and (14). 


Investigation of the Search Process 








The quality of the search process inthe i-th instant of time can be characterized by the search error Ay (i), 
equal to the difference between the value of y(i) and the minimum possible value of y at the relevant time: 


Ay (i) = u() — inj ty). (15) 


We introduce the mean error Ay mean(i) for the j-thsearch cycle: 


4 NG | 
Ay ()=sy > Ay. (16) 
i=3Nj 


mean 


From (10), we have 
N-1 


ay 0) =, Silla + Lj) + 3ka® + (Lj) + (3k + 1) a? 
k=0 
+[—a+L(j) + Bk + 2) a)*}. 


Hence, after some simple transformations, we obtain 





dy () = A[L3G) + Z(N*— 1a? + 3 (@—ay], an 
mean 
where L, (j) is the mean value of L(i) for the j-th cycle: 
Ly (j) =L(j) + “Sa. (1a) 


It is evident that Ly (j), and consequently Ay nean‘i), depends on the initial position L(0) of the system, and 
on the interference acting on the previous search cycles. We will obtain this dependence in explicit form. 
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From (5a), (13a), and (14a), we obtain 


Rlo(/) +5 + F(EM—aAW) 
OF+tM +n FEU)’ (19) 





2(j+1)= 
where L, (j) is given by the relation (18), and 


R= —2Aa(a—ay), 
S = aB, 
Q =2A (a —a)?, 
N-1 


EV) = + Soul, b), 
k=0 


N—1 


ni) = ar Y, bws(i—h, 
k=0 
2 N—1 
HC) = —> >} Sur (/, *). (20) 


k=0 
From the previously-indicated properties of the random interference § w(i), we have 


M {& (/)} = M {ny (/)} = M {p(7)} = 0, 
MEV) pV} = MEV 10)} = M (fn) p()} =9, 


MEU) =Min()} = +i MwD= >, (21) 


where D = M {8w%(i)}. 
Thus, 
MIE) +H) +10 = F- 


We will assume that for given A, a, N, and D, the quantity a is chosen so that 
D 
G> o (22) 
Then, with a probability close to one, the inequality 


QS lE() +n) +2)! (23) 
is valid. 


If we use the fact that the functions of the random arguments z(j + 1), that are the output of a real computer, 
are always bounded in absolute value, then the denominator of the expression in (19) can usually be linearized by 
the method in [7] in the neighborhood of the point —(j) = p (/) = n(j) = 0, and this yields 





eV +1) = (RL) +5+F 18) + nin [GF — POA + HO), - 


As will be shown below, the size of the test step a in an efficient system must actually satisfy the relation 
(22), and thus the investigation of processes in the system in question must be carried out on the basis of the approxi- 
mate equality (24). 
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When we substitute (24) in (12) and use (18), we obtain the recurrence relation 


Lj +1) = EL, (7) + C+ L,0)V()+W (), (25) 


where 
R+Q _ a 
oa et see 
C = SH3NQa _ aB+6NA (a—a)*a 
“3 Q ry 2A (a —a)* . 


Vj =-—GZRMO+a)+HM, 

W (7) = Wi (/) + Wa), 

Wi) = — Art ew — At hw — ee. 

Wa) = —aor (P+ PO) +EDHO—aMe OI. (26) 





We note the following properties of the random functions W (j) and V (j) that follow from (21). 
1. M{V(j)} = M {(W())} = 0. (27) 
2. Using the notation 

Dy = M{V*(j)},_ Dw = M {W*(j)}, Dyvw = M {W (j) VV ()}, 


we have 
Dy = For (28) 
Dw = M {W3(j)} + M (W3())} = (Be a il ay + ant (29) 
Dyw = M {W, (j)V (j)} = SBSP (30) 


QN 


3. From the above properties of the interference § w(i), for any n (n # 0) we have 
MV G)V(j +n)} = MW ()W( +n)} = MV GW +n) =0. (31) 


4. When we take into account the fact that L»(j) depends only on the interference of the previous search 
cycles, and does not depend on the interference of the j-thcycle, we obtain from (31) 


M {L, (j) V (j)} = M {L, (7) W (7)} = 0. (32) 
For j=1, let Lg(1)= Z%. Then, using the relation (25), we obtain 


L,(2)=C +W (1) + ZIE+V (pI. 
L,(3) = C {4+ [E+ V(2))} + W (2) + W (1) LE + V (2) 
+£1E +V (2) (E+V (4)I, 
L, (4) = C {4 + (E+ V (3)) + 1E +V (3) (E + V (2))} + W ()- 
+ W (2) (E+ V (3)] + W (1) (E+ V (3) LE + V (2)) 
+ ZIE+V (3) (E+V (2) (E+ (1). 


In the general case 
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Ly ())=CU+(E+VG—1)+(£4+V(—t(£+V(G—2))4+... 
+(E+V(j—11(E+V(j—2)]...1£+V (2))} 
+ W (j —1)+ W (j—2)(£+V (j —4)) 
+W(j—3)(E+V(—11(E+V(;—2))4+... 
+W (A)(E4+V(—1)1£+V (j—2))...[£+V (2) 
+ £(E4+V(j—1) (E4+V(j—2)])... 12 +V (AI. (33) 


By using the relation (33), we now find M {L,(j)} and M {L3(j)}. 
Taking ( 7) and (31) into account, we have 


j—2 
M {L, (j)} =C 5 E' + 2E*. (34) 
i=0 
If we now square and average the right and left-hand sides of (33), we obtain, after some transformations, the 
relation 
. j—2 ji-3 j-3—k 
M {Ly (j)} = C* >) ef + 2C (EC + Dyw) Sie* > E 
i=o k=0 i=0 
j—2 j-3 
+ Dw SE + £[ Lei + 2(BC + Dyw) 3 Ee], (35) 
i=0 i=0 
where 
e = E? ot Dy. (36) 


As can easily be seen, the necessary and sufficient condition for the convergence of the right-hand side of 
(35) for j > a is the inequality 


e< i. (37) 
Actually, when (37) is satisfied, the first three terms of the right-hand side of (35) are majorized by absolutely con- 
vergent series, and the term containing the factor £ , tends, for j > a, to zero no less rapidly than je. Thus, by 
using (37), we see that there exist finite quantities L, and L?, 
L, = lim M {L,(j)}, Lp = lim M {14(j)}, (38) 


independent of J and Z, and characterizing the control process in a steady regime. 
Substituting in (36) the values of the corresponding terms, we obtain 


a? 3 a*D 





diva aap + 2 Ate —ayn* (36a) 
It is obvious that the fulfilment of (37) is obligatory for an efficient system. Moreover, in order to obtain a 
su tapidly reversing search process, the quantity « must be chosen sufficiently small in comparison with 
unity. 


Let the range of possible values of « and A for a given object of control be given by the inequalities 


—%QAaQaty, = OC AS AS Amax (39) 
characterizing the control process in a steady regime, 
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action and the maximum possible value of € is € max We obtain the equation for the determination of the minimum 
possible value of a: 





Then, if we assume that the quantity ¢= is given, and starting from the requirements of rapid 


pe m ? “a (40) 
"(ra Ga 


“min 
In Fig. 2, we show a family of curves for (40), which can be used for the graphical determination of a,,4,, for 
given y, a», and €,.,,. It is obvious that the quantities A,j, and ain determine the minimum value of Q 
in (20). The curves in Fig. 3 were obtained from (40), 





and are a family of curves for Ss - as a function of 64/n 
gand €..,,- As can be easily seen from Fig. 3, the @ min r 
aur 7 


relation (22) must actually always be satisfied to obtain 
sufficiently small values of €,,,,, and moreover, in the 
a3t 
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region of small values of (gg < 1), the relation (22) is necessary for any €p,y (€max < 1). It should be noted 
here that the region of small g, as will be shown below, is more important in the application of the extrapolational 
search method. 


If we assume that (37) is satisfied, then we obtain the expression (38) in explicit form. If we use the fact that 
E? < ¢, then from (34) and (37) we obtain 


l= 5: (41) 


If we square and average both sides of (25), and use (32) and the fact that for j -- a, we have 
M {L,(j)} = M {Lp (j + 1)} = Ly, 
M {ZL} (j)} = M {Lp (j + 1)} =Z,, 
then we obtain 
[3 = F*L3 + 2CEL, + L3 Dy + 2L, Dvw + Dw. (42) 
Hence, taking (14) into account, we have 


v2 (1+ B)C*+ 2Dyy C + (1 — B) Dy 
Ly = (i—E) (i—E*—D,) , (43) 
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The substitution in (43) of the corresponding quantities in (20), (26), (28) - (30) yields 








6a*BD D + 3NB* 
ne Fags (0B + ONApAa)* (2p — a) + SEE, ii Natgt “Watg) we 
4a? [2 (2p — \— nal 


where p=a — a. 


From (44) and (17) we obtain an explicit expression for the mathematical expectation of the mean error Ye, 
over a cycle in the steady regime: 


Yex = limM {Ay (/)} = A[T3-+3(N?— thot +3 p']. (45) 


mean 


The value of Y,, obtained in this way can be used as a measure of the operational precision of a system. 


A Comparison of the Extrapolation and Gradient Methods of Extremum Search 





We now consider the object of control, described above, when the gradient method of extremum search is 
used. 


In this case, the control signal at the object input is given by the Eq, (8), in which x, (i) is given by 


, b for i= 2k, 
(=| to i=2k+1, = 


where b is the test step and k is the number of the search period (k=0, 1, 2,... ). 


The necessary changes in the signal x,, (i) are produced, in the general case, after each N search periods, i.e., 
at the times 2Nj, where j is the number of search cycles (j =1, 2, 3, . . . ). 


The control algorithm is given by the equations 





ty (j + 1) = zy (/) + 2 (7 + 4), (11a) 
S 
[wi (7, k) — wa (/, k)) 
Z(j+1)—*=2 aT ; (47) 


where B (B > 0) is some previously chosen constant, w,(j, k) = w(2Nj + 2k) and we (j, k) = w(2Nj + 2k + 1) (k =0, 
1, eee »N-1). 


The investigation of such a system was carried out by A. A. Fel*dbaum in[5]. We denote by Y,, the mathe- 
matical expectation of the mean of the tracking error over a search cycle in the steady regime. Then for the case 
N =1 and c,= 6 =0, we obtain, from the results of [5] and (8), 





: | : . Bbq? D 
Yor rand {= ty (27) + y (27 + 1)i} =A | 4a + RAa (Rh — Ag) a | , “ 


ee, (Al 4 
q=b 7 
Here the condition for convergence of the process will have the form 


i1— Stl <4, (49) 
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which corresponds to the system of inequalities 
gq=b—5>0, Aq< Bb. (50) 


Let the range of possible values of « and A be given by (39). For simplification of the description, we will 
assume further that the quantities a» and A;pin in (49) are equal to one; this assumption can always be made, since 
a transformation of the scales of x and y will always yield this result. 


In order that (49) be satisfied for the whole range of values of a and A, the quantity B, as can be seen from 
(39) and (50), must satisfy 


B= KAgax (1+), (51) 


where K(K > 1) is some given number determining the stability and rapidity of search for the extremum with A 
= Amax- 
Substituting (51) in (48), we obtain 


i 
2 on 
KAP a y( aa 55) a? 








Y,, =A te D 4+}. (52) 
gr a(14+5) swe KAmax(1 +35) -1) 
A(1 +35) 


It is apparent from (52) that, for A = A,,j,, =1 and ao =1, the inequality 


¥ i A? (53) 


gr max 
holds, where Yor max 48 the maximum possible value of Y 5, for a given range of variation of a and A, 


In order to compare (52) with the results obtained in the extrapolation case, we set N=1 and 6 =0 in (44) 
and (45), and obtain 


Da® D\} 
2a? (2p — 3 735 
a [7202 (2p — a) + ara + aca) | s2ph. ts 


[ep] 








in (45a) there is only one constructive parameter a, which is chosen by taking (40) into account for a given 
€max- As soon as the value of a is chosen, the expression in the curly brackets on the right-hand side of (45a) is 
a monotonically-decreasing function of A, and the maximum possible value of Y¢, is limited by the inequality 


Yex < AmaxK1, (54) 


where Kj is a constant, independent of A, equal to the maximum value of Y,, for A=1. 


It is apparent, from a comparison of (53) and (54), that for any values of the constants K and K, and for 
sufficiently large A,,,,, the quantity Y,, ,,,, cam be as much greater than Yex max a8 we please. We have thus 
reached the following important conclusion. 


For any level of stationary, random noise at the output of the object of extremal control being considered, 
and for a rate of drift a of the characteristic of the object in the range (39), where a, is any nonzero quantity, it 
is always possible to find a range of possible variations in the coefficients of the characteristic A of the object, 
such that the maximum, mean, steady error Y gr max in the gradient system will be as large as we please relative 
to the analogous quantity Y., m4, obtained for the extrapolation system, and in this case the rate of search in the 
extrapolation system, for all values of A and a, will not be less than the initially given arbitrary value. 
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It is evident that the smaller the value of K,, then the smaller will be that value of Ajax we may begin with 
for the extrapolation search to give a better result than the gradient search. 
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It is obvious from the expression (45a) that Ky decreases with increasing p, which, in its turn, depends on the 
size of the test step a. As can be seen from Fig. 2, a comparatively small value of a is obtained for the region 
0 < ysl, and for this reason the region indicated is the most suitable for applying the method of extrapolation 
search. 


Example. In Fig. 4 and Fig. 5, we show the dependence of Ygr and Yex on the values of A and a, obtained 
by calculations for the following conditions. 


For simplification of the calculations, it is assumed that the dispersion D of the random interference is propor- 
tional to A?, andso 6 = = =const. Two variants were considered: 6 = 0 (the curves in Fig. 4), and 6 
= 20 (the curves in Fig. 5). In both variants, the range of possible values of A and a were 


1<A<2, —1i<ac<i. 


The parameter B of the gradient system was chosen according to (51) for K=1.1. The size b of the test step 
in the gradient system was chosen from the condition that the minimum value of Y 9; max be obtained, and this 
yielded b = 5.4 and b = 5.3 for § =0 and 6 = 20, respectively. 


The size of the test step in the extrapolation system was chosen from the condition € max = 0.55, which, from 
Fig. 2 for a =1, gives a = 2.3 and a = 4 for 6 =0 and 6 = 20, respectively. 


The curves in Figs. 4 and 5 were obtained by substituting the above values for a = — 1, 0, +1 in the Eqs. (52) 
and (45a). 


The symbols “g* and “e" in the diagrams denote curves for the gradient and extrapolation systems respectively. 
The values a =+1, — 1, 0, correspond to the continuous, dashed, and dot-dashed curves, respectively. 


It can be seen from the curves that, in the case § = 0, the extrapolation system is more accurate than the 
gradient system for all values of A and a lying in the given range. For § = 20. the quality of operation of the 
extrapolation system becomes less good, and values of A and a can be found for which Yor > Ye,- However, as 
before, we have Y gr max > Yes max 


It is obvious that in both the cases we have considered the application of the extrapolation system ensures 
not only a more rapid search, but also gives greater accuracy in the steady regime. 


SUMMARY 
It follows from what has been said above that the benefit of applying the extrapolation search method depends 


on the size of the two ratios A,,4,/Amin and te a*/D. For sufficiently large values of these two ratios, the 
extrapolation search method yields not only a more rapid rate of search, but also a greater accuracy in the steady 
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h regime than the gradient method. From physical considerations, we must expect that criteria, similar to those in- 
troduced above, exist for the comparison of the gradient and extrapolation methods of search for cases when the 
output coordinate of the object is a second-degree polynomial function of several variables. 
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THE APPLICATION OF THE METHODS OF STATISTICAL DYNAMICS 
TO THE COMPUTATION OF THE CHARACTERISTICS FOR CERTAIN 
AUTOMATION OBJECTS 


S. Ya. Raevskii and N. S. Raibman 
(Moscow) 


Translated from Avtomatika i Telemekhanika, Vol. 22, No. 11, 
pp. 1466-1474, November, 1961 
Original article submitted January 20, 1961 


The paper studies the general problem in the theory of determining the operator for automation 
objects from the statistical characteristics of the input and output. The proposed theory is con- 
cretized for one class of technological process. An example is given involving the computation 
of the characteristics of an automatic assembly line. 


The methods of statistical dynamics which have been devised for automatic control systems have acquired 
wide application in other fields of science and engineering [1]. This can be explained both by the presence of an 
extensive class of phenomena in scientific and engineering problems which are statistical in nature and by the 
universality and generality inherent in the methods of statistical dynamics. 


One of the important problems which can be solved by the methods of statistical dynamics is the problem of 
determining the characteristics of objects with statistical inputs. Here the object proper can introduce statistical 
distortions due to certain internal sources into the input signal which it converts. The determination of the character- 
istics (in particular, the operators) for such objects is required in a great number of cases. A partial survey of the 
papers devoted to this subject is cited in [2]. 


In this paper we study general questions involving the theory of determining the characteristics of objects; 
then this theory is concretized for one class of technological processes. In conclusion we study an example in which 
the results of the theory are applied to the computation of one of the characteristics of an automatic assembly line 
which manufactures rings for conical ball bearings. A knowledge of this kind of characteristic is required both for 
increasing the operating efficiency of automatic ball bearing lines and for the over-all automatization of production. 


1. The General Theory fer Determining Object Operators from the Statistical 





Input and Output Characteristics 





Assume that the random signal X(t) is applied to the input of a certain object as shown in Fig. 1. As a result 
of this input a random process Y (t) is produced at the output of the object. The random processes X (t) and Y (t) can 
be measured. From the results of these measurements it is necessary to determine the characteristics for the objects. 
The characteristic of the object shall be defined as the operator by means of which this object is described. The 
result of applying the signal X(t) to an object with a certain operator can evidently be written in the form of the 
following equation: 


A,X (8) =¥ (2), | (1.1) 


where A, is the object operator and the subscript t indicates that this operator depends on the variable t. In the 
most general case there is no information on the operator which describes the object. In view of this, it is natural 
to assume that the operator A is random. Therefore, it is not possible to determine the operator proper, and we can 
only provide an estimate of it A,” which will then serve as the characteristic of the true operator A;. It is natural 
to require that the estimate of the operator be close to its true value in the sense of a certain criterion. This is 
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equivalent to requiring that the random function 

Y* (t) = A; X (s) (1.2) 
must be close to the random function ¥ (t). We know [1] that similar problems are solved in the following manner 
in statistical dynamics. We formulate a certain function p[Y(t),/’*(t)], which depends on Y (t) and Y* (t), and 
impose the requirement of a minimum on the mathematical expectation of that function; i.e., 


MplY(t), Y*(é)] = min. (1.3) 





The choice of the function p [Y, Y* ] depends on the criterion adopted. 
In order for relationship (1.3) to be satisfied it is sufficient that 


X (0 + | Object |> ¥ (oD 
Fig. 1. Object with a random 
input and a random output. MplY (2), ¥*(¢)/x(t)] = min, «.4) 











i.e., the function p [¥(t), Y* (t)) must be minimal for a fixed realization x(t) of the random function X(t), The 
condition governing the minimum of the relationship (1.4) is 


Fi 
ey” 





Mpe\Y (t), Y*(t)/2(t)] =0. (1.5) 


In the particular case where p [Y, ¥Y*] =(Y— Y* * (this corresponds to the criterion for the minimum mean- 
square error) we obtain from (1.5)* 


Y* (t) = A;X (s) = MY (t)/z (0). (1.6) 
Equation (1.6) makes it possible to determine the optimal prediction of the random operator A from the 
criterion governing the minimum mean-square error in the class of all possible operators. If we limit ourselves to 


the class of linear operators, it follows from (1.6) that it is not difficult to obtain another equation. In fact, we shall 
reject the idea of fixing the realization x(t) in (1.6). Then this relationship will be written as 


A;X (s) = MY (t)/X (2)}. (1.7) 
Multiplying both sides of (1.7) by X(u) and taking the mathematical expectation of both sides, we obtain 


M {A;X (u) X (s)} = M [¥ (t) X (u)], (1.8) 


since 


M {X (u) M LY (t)/X (u)l) = M IY ()X(@). 


In view of the linearity of A} the operator M is commutative with A‘ for the most general assumptions, 
Without limiting the general nature of our analysis it is also possible to assume M([X (t)}=0, M {¥ (D} =0, Then 
(1.8) is written as 

At K xx (u, s) = K xy (t, u), (1.9) 
where K,.(u, s) is the correlation function for the random process X(t), and Kyy (t, u) is the mutual correlation 


function for X(u) and Y(t). Equation (1.9) determines the optimal mean-square prediction of the random operator 
A, in the class of linear operators. A particular case of this equation is the integral Wiener-Hopf equation. 


* See also [1]. 
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Equations (1.6) and (1.9) determine the predictions of the random operator. It remains for us to clarify the 
problem involving the nature and quantitiative value of the error which is allowed when the random operator Ar is 
replaced by its nonrandom prediction At. For this purpose, we shall formulate the expression for the indicated error: 


e (t) = A, X (s) — A; X (s). (1.10) 


It is easy to show that the total error at the output of the object with the random operator A, is predicted on 
the basis of the dispersion D [Y (t)] and is written as 


D{Y (t)}) = D(Y*()) + Died), (1.11) 


where Y* (t) is determined in accordance with (1.6) and ¢ (t) is determined in accordance with (1.10). 


In fact, 
D{Y (t)] = M [A,X (s)]* = M {[ A, X(s) — A; X (s)] + Ar X (s)}? 
= M {e(t) + ¥°(t)}? = D[e(t)] + DIY" (t)] + 2M {e (t) Y° 0}. 


We shall now turn to the last term in the chain of equations cited above. In accordance with the chosen 
minimum mean-square error criterion we have the condition 


x 


aye MAY (t) —¥°(W)? = 0, 


M{Y ()] = M[Y* (i), 


L.e., 
M {A X (s)} = M {A; X (s)}. (1.12) 


Note that the latter relationship expresses the condition governing the unshifted nature of the prediction A¢ 
of the random operator A,;, since it is equivalent to M [e (t)] =0. Further, it is possible to write 


M {e (t) ¥* (t)} = M {A,X (s) A, X (s) — A; X (s) A; X (s)} 
= M {M [A,X (s) A; X (s) — At X (s) A; X (s)}/X (s)}. 


Since the internal operation of mathematical expectation does not affect X (s) it follows that due to (1.12) the 
mathematical expectation of the square brackets (and therefore of M[« (t) Y* (t)]) goes tozero, Thus, (1.11) has 
been proven. The physical meaning of this relationship resides in the fact that the over-all error at the output of 
the object consists of a) the error due to the random input signals corresponding to the fixed value of the random 
operator, and b) the error produced by replacing the random operator by this fixed value (as we have already noted, 
the fixed vakie of the random operator is assumed to be its unshifted prediction). 


2. The Application of the General Theory to the Prediction of Object Operators 





and To the Computation of the Errors of Technological Processes 





A typical object for a number of processes is the technological chain of sequential operations shown in Fig. 2 
(the process of machine design, instrument construction, metallurgy, etc.). In particular, such a technological chain 
may be equivalent to a certain automatic assembly line. Assume that the input of a chain consisting of n sections 
is subjected to a certainrandom factor X, (t) (our subsequent reasoning is valid when any infinite number of random 
factors operates). It is required to determine the effect of this random factor and the random signals produced by 
the sections proper on the spread of the output parameter X(t). These effects can be expediently considered in 
the following manner. 
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First, we shall study the effect of the input factor X, (t) and the entire technological chain on the output 
parameter X(t). In accordance with (1.10) and (1.11), we have 


Xn (t) AyXq (s) An X (s) -+- e) (t), 
D{X,,(t)| — D{Ay X (s)1 4 Die, 


where A,, is the unknown random operator for the entire chain of n sections, and A,,"is its prediction. 


A 
i oe 
—————as ee = 


= 
Ale) pelt aul) are 
Boe eee ee Oe Oe | 


Fig. 2. Block diagram of a production object. Fig. 3. Subdivision of the technological chain 
by means of the operators (k=1,..., 1m). 























Then we shall take into account the direct effect on X, (t) of the random factor X, (t) which acts at the input 
of the second section and the remaining portion of the chain consisting of the sections from the second to the n-th, 
We have 


X(t) = Aty Xj (8) = Ate X1 (8) + &2 (2), 
D(X, (t)] = D(A X1(s)1 + Dies (I. 


Continuing this process, as shown in Fig. 3, we evidently obtain the following two chains of equations: 


Xn(t) = At Xu-1(s) + ex (0), (2,1) 
D [Xn (t)] =D{[Aj),Xx- (s)] + D [ex (t)| (k = 1, 2, -o- BM, n), (2.2) 


Summing Eqs. (2.2) with respect to k, we obtain 


2 = > dD, [Xn (t)] = 2 dD [ Aj, Xn (s)] + 2 D [&% (t)}. (2.3) 
=] [== 


k=1 


The first of the sums on the right side of (2.3) takes into account the effect of the errors at the inputs X, (t) on 
the output parameter X,(t); the second sum takes into account the resultant effect of the random perturbations in 
the technological chain proper. 


We shall now proceed to a computation of these sums, As practical experience shows, for specific computations 
of a number of processes which arise in a number of different engineering fields we can make the following assump- 
tions; 


1) All the errors X;, (t) (k=1, 2,..., 0-1, n) are delta-correlated; i.e., 
K xx (ty, te) = Smid (ty — te); (2.4) 


2) All the sections of the technological chain are linear; ie., 
T 
AtmXm (8) = \ W(t, s) Xm (8) ds. (2.5) 
0 


Note that writing the operators for machine-design objects in the form (2,5) does not contradict the discrete 
nature of machine-design processes, since this is easily taken into account using the special form that involves writing 
W(t, s) in terms of a delta-function [1]. However, this form of notation will not be required for our subsequent 
analysis, 
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We shall now determine the unshifted predictions A*,, of the operators Ar. 


In view of the linearity of A, we have the eg result in accordance with (1.9), (2,4), and (2.5): 


T 
( W's, (t, s)p(u—s) ods = Kt, u). (2,6) 
i) 


From this we obtain 
e i 
W(t, ¥) == KS) (t, u) (2.7) 
Sin 
by using the known property of a de’:a- function. 
This determines the prediction of the operator, and we can write 
T 


At X(s) = | K% (t, s) X (s) ds. 


mo 


(2.8) 


We shall now compute the first of the sums on the right side of (2.3): 


= 5} Dl AuXr— (s)] = 


k=1 


0 
TT 
x M[X (s) X (s,)] dsds, = 5} 7 \ Na 
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X p (8, — $2) dsyds, = 
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In particular, for constant mutual correlation coefficients we have 


n 


D, = >} 


k=1 %%—1 


[Koy ’P. 





(2.9) 


In order to compute the second sum we note that in view of (1.11) and (2.8), we have 


D{ex] wy (4m _— ko? — > ) Xin (s) | 


= Dr{Xn(s)]—D [ Rav" ao Xn *) |= Dy {Xn (s)] — (Ky PG 


o2 

Sh 

2 

k—1 n-k— 


From this we obtain 


D (en) = Dr (Xn(s)]—ontrSy P, (2.10) 
where - » is the coefficient of correlation between the output and the k-th input. 


Making use of (2.10), we can Poe obtain 


i 
| 
Si 
Mi 
, 


Dy = ¥) DalXel— west ot) = of (1S es" P). (2a 


k=1 
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The use of statistics dynamics methods makes it possible to determine the values of Dyn at the output of an 
automatic assembly line from the values of D,, for each of the processes used along the line. Thus, for the assump- 
tions made above, D, yn can be determined from the formula (3, 4): 


n n—1 n—1 
Din _— > Aik + 2 b> >; Ai. AuK gx, (2,12) 
k=1 k=1 l=k+1 
or 
n 
Din _= >> Auky,x,, ’ (2.13) 


i=1 


where A‘, are the predictions of the operators for each of the processes used along the automatic line (Fig. 4). 
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Fig. 4. Scheme for processing the vibrating track 

for the outer ring of a *7815K1" ball bearing on the 
automatic line of the First State Bearing Plant. 1) 
Lathe processing (8-spindle double-block automatic 
lathe "02S01"); Il) heat treatment @OKB-134" furnace, ; 
Ill) preliminary grinding (automatic internal grinder 
*02S533B"); IV) finishing grinding (internal automatic 
grinder "02S33"); V) fine grinding (internal automatic 
grinder “02S33A"). 


3. An Example Illustrating the Computation of the Characteristics of an 





s.utomatic Production Line 





As an example we shall study the problem of determining the characteristics of the automatic production line 
for the First State Bearing Factory; the line processes the outside ring of the *7815K1" ball bearing. We shall study 
the processing of one of the basic parts — the vibrating track.* The scheme showed the processing of this part 
on the automatic production line is shown in Fig. 4. 


The manipulation of the experimental data consisted of predicting the values of mathematical expectations, 
dispersions, and correlation functions for the realization of each of the processes: it also included the determination 
of the mutual correlation functions and the coefficients of correlation between the individual processes performed 
along the automatic line in order to produce the vibrating track. 


Empirically normalized correlation functions, as computed from the indicated experimental data, are cited 
in Fig. 5 for the results of the operations involving rolling (preparation), heat treatment, and grinding. The values 
of the empirically normalized correlation functions can be taken as approximate predictions of the normalized 
correlation functions for the processes under study. From Fig. 5, it is evident that for the very first shift by one 
part (n,; = 1) the prediction of a correlation coefficient is close to zero, and therefore the errors for the investigated 
batches are delta-correlated. 


When the volume of investigated parts in a batch is considerable (i.e., for appreciable realization length) the 
determination of the basic characteristics of the processes and the interrelationships between them requires an 
appreciable amount of computation. Moreover, the memory of modern electronic digital computers is inadequate 
for large realizations. A substantial reduction in the volume of computations can be achieved by presenting the 


* The characteristics of the line are computed from the results of measurements performed by the Precision Laboratory 
of the Machinery Institute of the Academy of Sciences, USSR. The authors express their appreciation to M, 1, Kochenov 
and V. I. Sergeev for supplying the original materials. 
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data in correlation tables for which methods of formulation have been well-developed in textbooks on mathematical 
statistics (5-7). 


The formulation s cnmaienen tables for random functions of stationary processes is performed on the basis 
s the counseed values x() and x (t +m) for the computation of correlation functions, and according to the values 


x (t) and Y (t + m) for computing the mutual correlation functions (m =0, 1, »k sw n/4, where n is the 
volume of the batch). The number of correlation tables depends on the number m (i.e., it depends on the length 
of the realization). The complete processing of the correlation tables is performed on an electronic digital computer. 


As’an illustration, Fig. 6 shows the correlation field, 
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operations: a) rolling; b) heat heat treatment (y); M{y/xz} = 1,034% — 3,495; 
treatment; c) finishing grinding. n = 0,751; r= 0,750 


In view of the fact that the errors are delta-correlated, it is possible to limit ourselves to a correlation field for 
just a case where m = 0, in the example under study. 


From each of the correlation gables we determine the values of the normalized mutual correlation functions 











K,.,, (t) 
in xy 
Ixy (t) V D, © D, © (3.1) 
and the values 
M{MIY —_— 
n()=Y vor (3.2) 
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where n (t) is a natural generalization of the concept of the correlation ratio for a stationary random process which 
is treated as a time function. For delta-correlated cases, we can limit ourselves to determining r(t) and n(t) for 
m = 0 in the example under study. 


From the values of the correlation coefficients r and the correlation ratios n computed for each value of m 
it is possible to establish the linearity of the investigated process by determining the criterion [5-7]. 


a =) Ale —*) 
P= 3.3 
G—Dd-w) = 
where n is the number of experiments which are performed; h is the number of intervals of the correlation table 
which is used to perform the computation of the values of n and tr. 





Fox each case, the value of F computed according to (3.3) is compared with the theoretical value F_ for numben 
of degrees of freedom equal to (n - 2) and (n - h). If F < Fy then the hypothesis concerning the linearity of the 
system does not contradict practical experience and can be accepted. The hypothesis concerning system linearity 
can be verified for the functions n(t) and r(t) on the basis of the maximum difference between these functions for a 
fixed argument t. The values of Fy are tabulated for definite values of reliable probability. For the example under 
study, the interrelationships between the heat treatment and lathe processes are the following: n=188,h=7, n = 
=0.751 and r=0.750. 


Then, in accordance with formula (3.3) F is equal to 6.0 if we exchange the numerator and denominator. For 
a reliable probability of 0.99 we have the value Fy = 9.08 [5-7]. 


Therefore, the hypothesis concerning linearity of the heat treatment process does not contradict practical 
results and can be accepted. We verify each of the operations in analogous fashion. 


The basic characteristics of the processes involved in producing the vibrating track on the automatic production 
line are shown in Tables 1 and 2 for two batches of parts on the basis of the indicated experimental results. 


From the specified correlation coefficients (Tables 1 and 2) we obtain the following values of D, and D, (2.9) 
and (2.11) for the example under study: D,p = 2160-10-? mm*, Dp, = 540-107? mm* 


Analogous results can be obtained for Dy, from formulas (2.12) or (2,13); for this purpose it is first necessary 
to determine the values of the predictions of the random operators A;}, for each process (see [3] and Fig. 4). Asa 
result we obtain: 


Aj, = 23-10-4, Aj, = — 306-10-4, Aj, = — 250-10-4, 


Aj, = 6.5-10-6, Aj, = 8690-10-*. 


The resulting operator predictions A*,, can be used to compute the output accuracy characteristic for the 
second continuous batch. In that case the values of the optimal operators A,), are associated with the characteristics 
of the second continuous batch, As a result we obtain og, = 13.67 p at the output. This value is in good agreement 
with the actual value of o, for the second continuous batch whose 95 % reliable interval is 12.4 < og, < 15.2. 


SUMMARY 


The application of the general methods of statistical dynamics to the computation of the characteristics for 
technological processes was illustrated in this paper on the basis of linear processes whose eriors are delta- correlated. 
For more complicated cases, it is possible to use the general devices and formulas cited in section 1. The paper 
examined only processes with one input and one output. It is not difficult to show that the methods cited above can 
easily be generalized to apply to multidimensional cases where 1) we are dealing with processes having many in- 
puts and outputs, and 2) the processes proper consist of several controlled systems. 
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THE SEARCH FOR EQUATIONS WHICH DETERMINE THE RELATIONS 


EXISTING WITHIN COMPLEX OBJECTS 


V. P. Borodyuk and G. K. Krug 


(Moscow) 


(Translated from Avtomatika i Telemekhanika, Vol. 22, No. 11, 
pp. 1475-1481, November, 1961 
Original article submitted March 3, 1961 


We examine the question of formalizing the control process in complex objects by means 
of the probability equations of statistics. We classify the various groups of parameters acting 
upon the object and examine the form of the statistical relations between the groups. We discuss 
the possibility of using the system of probability equations to increase the production efficiency 
of the objects. We present an example in which the process of working out the system equations 
by means of the method of correlation analysis, using the “Ural-1" digital computer, is described. 


1. Characteristics of the Complex Object 





There exist at the present time, in the various branches of industry, a large class of objects which are 
characterized by a multiplicity of interacting parameters. The physical nature of these parameters and their 
quantitative characteristics vary greatly and are partly unknown. 


We can divide the parameters acting upon such an object into a series of groups depending upon their character 
and the part which they play in the process. In Fig. 1, we have divided these parameters into four groups. The 
groups of parameters possess the following characteristics. 


4% 











sf K, 
; : 
ae object M3 
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Fig. 1. 


1, Group X parameters are input parameters. By input parameters we 
mean the totality of indices characterizing the quality of the initial products. 
The quantitative values of the input parameters are usually known: they are 
the data obtained by means of laboratory tests, chemical analyses, data ob- 
tained from various analyzing devices, etc. 


2. Group K parameters are output parameters. The output parameters 
contain information regarding the qualitative characteristics of the output 
product and also data regarding whether the object can be readily produced in 
the factory, the cost price of the product and other generalized indices which 
evaluate the accuracy and effectiveness of the object's operation. 


3. Group Y parameters are perturbation parameters. The presence of 
perturbation parameters vastly complicates and muddles the general picture 
of the process. The temporal characteristics of the perturbation parameters, 


their points of application to the object and the intensity of their action have a random character. Their material 
nature and qualitative meanings are unknown. In this group we include the uncontrollable indices of the properties 
of the input parameters, the action of chance (random) external factors, changes in the characteristics of technological 


equipment, etc. 


4. Group Z, the group of control parameters,characterizes the state of control actions, which are available 
to the operator. These include the indications of the cost-measuring devices, the positions of the control and re- 


gulatory devices and valves, etc. 


The model described may be applied to many technological processes in the chemical, metallurgical, paper 
and cement industries, Over and above this the model may be applied to the work of an entire enterprise (or in- 
dustry) where we can consider as the output parameters the profitability, unprofitable expenses or other generalized 
economic indices; we may in a similar manner choose the input and control parameters. 
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By varying the values of the control parameters the operator in charge of the process in the object tries to 
neutralize the perturbing actions of the uncontrollable factors and to achieve optimum operating conditions. 


The multiplicity of parameter interactions and the presence of many external perturbing factors greatly com- 
plicates the development of technological control programs and makes extremely difficult the design of programming 
cards, Therefore the operator must often rely upon his personal experience, and the effectiveness of his actions 
depends to a great degree upon the depth of his experience with the process, i.e., his knowledge of the process, and 
the sharpness of his intuition. 


2. Problems Involved in the Formalization of the Control Process 





The dominant role of subjective factors in determining the operation of a process prevents operation under 
optimum conditions for long intervals of time and leads to confusion, unstable operation and sharp variations in 
productivity; this is borne out by actual experience in dealing with complex objects. 


In order to increase the production efficiency and to obtain more orderly methods of control we must deter- 
mine the character of the relations between the groups of parameters and know the basic influence of the separate 
parameters upon the course of the process. We can obtain this information from a system of equations which for a 
given process or mode of operation determines the relationships between the various groups of parameters. 


It is obvious that in the absence of the uncontrollable perturbing forces Y, we could give equal to weight to 
the various functional bonds. 


The presence of intense random disturbances forces us to examine the problem from a statistical (probability) 
viewpoint. 


When describing the dynamics of the majority of complex objects subject to slowly changing, continuous 
processes the equations of the process may be written as a series of links associated with pure delays. In those cases 
where the delays are associated with transportation and shifting of the products of the process, determination of the 
time delay is usually not associated with great difficulty. In more complex objects the equivalent delay time 
between a pair of parameters can be determined approximately by means of the mutual correlation function ob- 
tained from the experimental data. In this process we take the delay time to be the average time for which the 
mutual correlation function computed for various intervals is a maximum, i.e., where the correlation between the 
pairs of parameters being studied is a maximum. 


In order to obtain the above-mentioned system of statistical equations for the relations between the parameters 
we must overcome significant difficulties inasmuch as the specifications for complex objects in the majority of cases 
make it much more difficult for us to use the method of artificial perturbations in obtaining the necessary character- 
istics and limit the application of physico-chemical modelling. The most effective method in the given case is the 
method of working over of the statistical data which is associated with the course of the process, in particular the 
method of correlation analysis. 


The method of correlation analysis, which analyzes a large number of experimental data, permits us to derive 
equations which show how the function varies on the average as one or more arguments of the function change while 
the other variables (arguments) remain unchanged. 


The success experienced in the application of the method of correlation analysis depends to a great degree 
upon the volume of statistical data available and the amount of information available within the statistical data 
regarding the problem in which the investigator is interested. As part of the method we can use the results of specially 
designed experiments and other technical material which has been documented and is stored in the archives of the 
industry. 


In order to clarify what we have outlined above we will examine the various variations possible in the relations 
between the various groups of parameters of the complex object (Fig. 1). 


1. The dependence of the group K parameters upon the input X parameters is determined by probability 
equations of the form: 
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+ DayjXF (t — tj) +... (1) 


J 


where s = ci, m=-n, moj, n<j, r= ci, us-vu-ew, u<j, v<j, w <j, K is the mean value of 
the i-th output parameter, X¢ is the input parameter, j is the number of input parameters included in the given 
equations, Ca is the number of combinations of j with respect tog and re is the delay time for a given channel. 


The form of the relation can be either linear or nonlinear, Usually, in practical examples, it is not necessary 
to consider components of higher than the third order. 


Such a system of equations enables us to predict the most probable values of a parameter K; in which we are 
interested, for given values of the input parameters and also enables us to evaluate a given situation. 


In addition, for a given amount of variation of the output parameters we can work out the limitations upon the 
allowable variations in the values of the input parameters. 


However, we cannot use these equations for the direct control of a process since the values of the input 
parameters are not subject to control. 


2. The dependence of the group K parameters upon the Z parameters is expressed by the equation 
K;(t) = \ize—u. (2) 


The probability equations of the system are analogous to Eqs. (1). 


The equations giving the dependence of the K parameters upon the Z parameters can be considered as auxiliary 
equations in the work of the operator. However, they do not reflect the role of the input parameters, which exert 
a basic influence upon the process. However, we may use this system in the following manner: If we substitute in 
the necessary equation the value of the parameter Kj corresponding to a given situation we can find the minimum 
value of the outlay for expensive products, i.e., we can determine for the time under consideration the most probable 
economic mode of operation. For the case where the results of the statistical equations has an extremal character, 
the optimum value of the control parameters may be determined from the null part of the derivatives corresponding 
to the probability equations. 


3. The dependence of the control parameters Z upon the X parameters is of the form 


. 


Z(t) = \ [X(t—1)].: 


The structure of probability equations for this system is analogous to that of the Eqs, (1). 


The data obtained from the equations permits us to find the most probable values of the control parameters 
for a given combination of input parameters. In each region to the values of the input, K, parameters there will 
correspond a system of equations connecting the X and Z parameters. The system of equations which gives us the 
desired values of the K parameters in a given region is of practical value to us. These equations permit us to pre- 
dict the values of the control parameters which will lead to the desired effect in a given concrete situation; the 
desired effect may be, for example, high productivity, good quality, economic production, etc. 
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We must bear in mind that in order to obtain the desired system of equations we may use that statistical 
material for which the desired effect takes place, for example,in the organization of the experiment or the investi- 
gation of available statistical materials (technically-documented data) we must choose only those values of input 
or controlled actions for which the values of the output parameters fall within indicated limits, 


4. The dependence of the group K parameters upon the simultaneous actions of the input and control parameters 
is described by probability equations of the form 


p=s 


K; = by + LbusX; (¢ —t;) + Ds OapX en (t — tm) Xn (t — Tn) 


p=1 


+ Y}bgpX3 (t — 13) +... +d, + DidyZ; (t — 1) 
j f 


q=l 
+ >} degZn(t — tr) Zz (=e) + Di dali (t — Tj) + noes (4) 
q=1 
where s = ci, m=sn, moj, n<j, l= ci, hy, h<f, gs, K is the average value of the i-th output 


parameter, j is the number of input parameters in the given equation and f is the number of control parameters 
in the given equation. 


Obviously, the given system of equations takes into account the action of all the control parameters, and gives 
a more complete mathematical description of the object being studied. Taking into account all the probable effects 
described in the cases above the given system of equations can be used to control the object for all values of the 
output parameters. 


Any probability equation has an empirical character and therefore it can be applied only over that region for 
which it was obtained (that is only over a certain range of values of the variables), Usually we are not justified and do 
not have a sufficient basis for extending its application beyond these limits. 


The larger the number of independent variables included in the equation the more accurate will be the control 
function and the smaller the deviation of the values obtained in the process from the values predicted by the 
equation. However, it is not expedient to include all the possible variables in one equation since then the equation 
becomes unwieldy and difficult to use. It is more convenient to divide the variables into separate subgroups, so that 
each subgroup is unified by parameters which are subject to a mutual influence but which are not subject to the in- 
fluence of the parameters in the other subgroups [1]. Such a division of the parameters simplifies the calculations, 
and leads to a more accurate model and presents a truer picture of the process being studied, 


Thus, we obtain a subsystem of equations which connects the given parameter K; with the various subgroups 
of controlled parameters. In determining the optimum mode of operation of the object with respect to the parameter 
Kj we must solve all the equations in the subsystem simultaneously. 


3. Mathematical Method of Determining the Equations which Express The 





Relations Between the Parameters 





The determination of the probability equation reduces to the determination of the coefficients of the equation, 
The method of correlation analysis makes use of the method of least squares. Of all possible values of the coefficients 
we use those which satisfy the condition 


p= X(Ky — K;)? = min, (5) 


i.e., the sum of the squares of the deviations of the experimentally-obtained values from the values calculated by 
means of the equation must be a minimum. In Eq. (5) ¢ is a function of the desired coefficients. We know that we 
can find the minimum function by equating to zero the partial derivative taken with respect to each variable. 
Going through this process for expression (5) we obtain m equations with m unknowns, where m is the number of 
coefficients in the probability equation. Simultaneous solution of these m equations permits us to determine the 
coefficients of the probability equation. 
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Having determined the statistical relation between the various parameters in the probability equation we must 
evaluate the degree or tightness of fit of the equation (2). This will permit us to determine the influence of the 
parameters which do not enter into the equation. In other words, the degree of fit of the equation determines the 
degree to which the equation holds in the presence of manifold disturbing influences. The tighter the fit, the more 
accurately can we predict the values of the function. 


The quality of the fit is expressed by the correlational ratio 
y= V* ‘ (6) 


Here o” is the dispersion which characterizes the total variation of the values of the given parameter Kj in 


yxy - 


the vicinity of its mean value K, — 8? is the dispersion of the particular means; it characterizes the 
n 





variation of the particular mean values of the parameter K, in the vicinity of the mean value K;. 


The physical meaning of the correlation ratio is described as follows: It ifdicates which portion of the total 
variation of the parameter Kj is due to changes in the arguments which are included in the given equation. Thus, 
in practice, n is a measure of the tightness of fit of the relation. 


The value of n lies within the limits 0 <n <1. If =1 then the relation is functional, i.e., all the parameters 
which depend upon the parameter K; have been accurately determined, 


If n =0 then there is no correlation between the parameters under study. 

If 0 <n <1 then we speak of weaker or stronger correlation. 

In the case where the correlation is linear the coefficient R may be used as a measure of the tightness of fit. 
The physical meaning of R is the same as that of n [2]. 


4. Example of the Determination of the Probability Equations 





The object studied will be a complex technological process, a continuous chemical process characterized by 
19 input parameters, five control and one output parameter, the output parameter being the quality of the final 


product. 


The statistical materials were gathered during the 
natural course of the process and they were worked over 
by means of the methods of correlation analysis. The 
control parameters were referred to the input and the 
process under study was considered to be an object with 
24 input parameters and one output parameter (Fig. 2). 
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Fig. 2. Fig. 3. 


All the input parameters were divided into subgroups (two parameters in each subgroup) in various combinations. 
This method of combination, as we have indicated in the article, permits us to take into account the correlation 
between the individual pairs of input parameters. We excluded from the possible variations those pairs which 
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included parameters which were independent or where possible dependence relations did not exist, We were thus 
enabled to reduce the total number of equations in the system to 66. We also assumed that the relation between 
each pair of input parameters which we were seeking was linear. In this process we bore in mind the following con- 
siderations. 


1. Detailed study of the technological process and of the results of physico-chemical investigations show that 
changes in the value of the output parameter as a result of changes in the input parameters are not critical. 


2. The range of changes of the input parameters comprises an insignificant portion (5-12 %) of the entire 


scale and we can expect that we can approximate the dependence in this portion by means of a straight line without 
incurring a large error. 


Thus, the desired system contains 66 equations of the form 
K = ao + 0X4 + aX, (7) 


In contrast to other similar works, for example the work described in [3], in the given case all necessary cal- 
culations were carried out ty means of the “Ural-1" digital computer. 


The preparation of the material for machine use (establishment of the scale, the making of the perforated 
tape, etc.) and the programming of the material took an intermediate classification programmer 5-6 working days. 
The machine time required for solution of the problem was 4.5 hours. 


The basic stages in the solution are illustrated in Fig. 3. 


Let us examine the actions performed in each stage: 1) obtaining the necessary initial statistical material; 
2) calculation of the means; 3) calculation of the dispersion; 4) determining the coefficients of the equation; 


5) calculation of the individual coefficients of correlation; 6) calculation of the coefficients of multiple correla- 
tion. 


The results obtained contain the following information: 

1) the coefficients of the 66 chosen equations; 

2) the coefficient of multiple correlation RKX;X j for each equation, These determine the tightness of fit of 
the relation between the given pair of input parameters X;X; and the output parameter K. 


In addition to the calculation of the partial coefficient of correlation RX Xj° which evaluates the significance 
of the proposed correlation between the given pair of input parameters, we must also determine the coefficients of 
correlation Ry, and Ryy,, which permit us to determine the influence of the associated input parameter upon the 


X 
quality of the eid product, We assume that all the input parameters involved in the correlation are independent of 
each other. 


The value of the coefficient of multiple correlation RKX{X; lies in the interval 0.01 - 0.34. 


After verifying the statistical significance of the coefficients of multiple correlation which we had obtained 
we chose six equations where the values of RKX{X; were within the interval 0.31 - 0.34. 


SUMMARY 


The existence of a system of statistical probability equations for a complex object permits us to determine the 
parameters which have a basic influence upon the course of a process and also the parameters of secondary importance 
whose influence is not of such significance or does not affect the process to such an extent, and also to evaluate the 
situation for any existing conditions and therefore to develop a system of tolerances with regard to the quality of the 
inputs. 


Using the system of probability equations we can make up cards for the technological operations which will 
greatly simplify the control of complex objects. These equations may be used as the basic equations in the con- 
struction and design of complex automatic systems which require multiple-relation stabilization. 
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Finally, if the system of equations which we have obtained describes the operation of the entire process it will 


help basically to evaluate the technological-economic indices of the operation and permit a rational distribution of 
the means and materials. 


1, 
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Theoretical basis is given for the use of the describing function in approximate determina- 


tion of periodic states in nonlinear pulse systems of various types of pulse modulation when the 
nonlinearities of the pulse elements are of arbitrary character. 


INTRODUCTION 





The describing-function approach introduced into automatic-control practice by L. S. Gol'dfarb (1-3) is here 
extended to nonlinear pulse systems. 


Nonlinear pulse automatic systems (NPAS) embrace a wide variety of devices from ordinary sampled-data 
control systems to those which include a digital computer (4, 5). Periodic motions may arise in an NPAS as well as 
in an ordinary, that is not a pulse nonlinear system. To determine fully such motions in the case of high-order 
systems proves impossible as no analytic methods are. available to solve nonlinear difference equations [5]. 


At first, simple NPAS were studied. A relay-pulse system without dead zone was in particular analyzed in 
[6] by means of difference equations, Periodic states, however, were not considered at all. 


Work[7]dealt with approximate determination of periodic states in relay- pulse systems; the pulse train pro- 
ceeding to the continuous part of the system was replaced by an equivalent sinusoidal area. 


The describing-function method was first applied to analyze an NPAS in[8). There,a relay-pulse system with 
dead zone was investigated. Exact methods were found to study periodic states in relay-pulse systems with no dead 
zones. These methods were developed by Yu. V. Dolgolenko [9] who used discrete Laplace transforms and by V. P. 
Kazakov [10] who used the frequency technique. The latter showed that by neglecting higher harmonics of the exact 
solution, one obtains results similar to those of Chow for relay- pulse systems without a dead zone. In[11] the de- 
scribing function approach was used to analyze periodic states in width- pulse systems. 


In the above investigations the describing function was used in various ways in NPAS and proved applicable in 
those special cases; however, no universal methods emerged as regards its application to nonlinear pulse systems. 


It is known that there was one attempt to systematize the use of the describing function for NPAS, Namely, 
R. E Kalman [5] in discussing different methods of analyzing a NPAS refers to the cumbersome exposition of (8); 
he comes to the conclusion that the obstacles are insuperable when trying to apply the describing function to some 
more involved cases, 


It is envisaged in the present paper to give an outline of the main features of the describing function in the 
analysis of a NPAS. General theoretical propositions developed by the author form the basis for the determination 
of periodic states in systems containing a digital computer.* 


* The subject of the present paper constitutes the theoretical part of “Periodic states in systems with a digital 
computer,” a contribution of the author delivered at the Sixth Conference of Young Research- Workers of the Institute 
of Automation and Remote Control of the USSR Academy of Sciences on January 19, 1959, It was published in a 
condensed form in [12], 
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1. General Comments 





The essential difference between the nonlinear pulse systems and the ordinary continuous nonlinear systems 
lies in the signal pulse modulation. The distinction manifests itself in that in addition to the usual nonlinearities of 
a control system such as saturation, dead zone, backlash, etc. an NPAS can also contain special "pulse" nonlinearities, 
in particular: the nonlinearities of the pulse elements. 


A NPAS shall be considered in the following such that it can be subdivided into a nonlinearity, a pulse one or an 
ordinary one, and into a linear part. 


Figure 1 represents a simple NPAS whose characteristics re- 
semble those of ordinary nonlinear systems, In this system the 
pulse modulation is implemented by the linear pulse element LPE 
which divides the linear part of the system into two linear filters 
LPt, and LPt,, The nonlinearity is represented by an ordinary non- 
linear element NE. 











The transformation of the signal by the linear pulse element 
brings about the occurrence of high-frequency components in the spectrum In this sense the LPE as well as the NE 
acts as a “generator® of higher harmonics.* The linear filters LPt, and LPt, suppress to some extent the high- 
frequency components generated by the NE and the LPE. Therefore, the application of the describing-function 
method to such systems can be based on assumptions similar to those made when determining the periodic states of 
ordinary nonlinear systems with several nonlinearities (3). By these assumptions both filters, LPt, and LPt,, are 
supposed to transmit only the fundamental harmonic of the periodic process taking place at the output of the NE 
or the LPE respectively. The equivalent complex-valued gain of the nonlinear element depends here not on the 
magnitude of the frequency of the fundamental harmonic itself but on the ratio of this frequency to the repetition 
frequency of the pulse element, that is,it depends on the relative frequency. 


In contrast to the previous case, in an NPAS with a pulse nonlinearity, there is no filtering of harmonics 
between the nonlinear and the pulse conversion,and all the high-frequency components generated by one conversion 
are submitted also to the other one. As a result of this there appears a signal at the output of the pulse nonlinearity 
in which the components of composite frequencies play an essential role. 


In the application of the describing-function approach to this type of systems, the composition of the higher 
harmonics from the nonlinear and pulse transformations is of great importance. 


Further on, we shall be extending the describing-function method to systems with pulse nonlinearities, or more 
precisely, to systems with nonlinear pulse elements. 


2. Statement of the Problem 





A nonlinear pulse automatic system (Fig. 2) with a nonlinear pulse element is now considered; we take f(t)=0 
(€=t/T, is the dimensionless time). The system is represented by a closed network consisting of a linear part LPt 
which can be described by a linear differential- difference equation, and of a nonlinear pulse element NPE. 


A train of modulated pulses z(t) (Fig. 3) proceeds to the input of the LPt; the pulses may be modulated in one 
or more ways (amplitude — width — position or phase, shape-modulation etc.), There is a nonlinear dependence of 
the pulses z(t) on the discrete values u[n] of the input signal u(t) to the pulse element. It is understood that this 
nonlinear dependence is such that any discrete value identifiable with its conventional “short” pulse is converted by 
the pulse element, as if it were a nonlinear pulse shaping device only, into a specified pulse at the output. The 
actual form of the nonlinear shaping device is quite arbitrary. 


The nonlinear pulse conversion takes place, as a rule, simultaneously with the variation (depending on the 
magnitude of the input signal) in the time-dependent parameters of the pulse element itself. Therefore, the usual 
representation of a pulse nonlinearity (by means of a static characteristic, such as a known time delay, etc.) is 
bound to meet with the difficulties, which are well-known, and it does not, generally speaking, make sense. A 


* The two generators, however, differ in the manner in which the spectrum is formed. 
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pulse nonlinearity, in contrast to the above, is characterized by the modulation characteristics on =F,, (u(n)) of 

the pulse element,which relate the discrete values at the input to the values of the pulse parameters at the output 
(see, for instance, Fig. 4 where the amplitude-modulating characteristic is shown; in this case on = ip, where ap 

is the dimensionless amplitude of the n-th pulse). In addition, a pulse nonlinearity is characterized by the repetition 
frequency w», = 21/T, of the pulse element as well as by the nominal amplitude Kp of the output pulse. 


There arise periodic states in the systems under consideration due to nonlinearity of the pulse element, The 
frequency w,= 2m /T, of the periodic state is such that the repetition frequency ws, is a multiple of the former. 
Thus, when trying to determine the “free* periodic motion of a NPAS, one encounters a nonautonomous oscillating 

system due to time quantization. Also the periodic motions 
in such systems are not self-oscillations in the usual mean- 
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Subsequently, we shall be interested only in NPAS such that the approximation due to the use of the describing 
function is valid; we shall determine approximately the periodic states in such systems. It will be assumed for the 
sake of simplicity that the periodic motion has no constant component. 


3. Existence of Periodic States 





In accordance with the above-made stipulations, a periodic state is considered such that 
pw AY, 
N= hg P= 2 (1) 


is an integer, and also such that at the input of the NPE the periodic state is approximately determined by its 
fundamental harmonic 


u = Ccos(+ t+ y) (_4<¥<7)) 


where @, = = = @,T, is the dimensionless frequency of the required periodic motion. 


We write the complex-valued amplitude of the input signal into the NPE as 


U = Cel, (2) 


and the complex-valued amplitude of the fundamental harmonic of the NPE response to the (2) signal we denote by 
Z,. In accordance with this notation the equivalent complex-valued gain of the NPE is now 


=H 


i 
According to the describing-function methods [1, 2, 3], the equation of the periodic state is approximated by 
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[1 ms Kir |v =0, (3) 


where K (j@) is the frequency characteristic of the system linear part and W = wT, is the dimensionless frequency. 
It follows from Eq. (3) that 


K(j 7) =— ; (4) 


is the approximate condition for the periodic state to exist. 





The approximate-condition formula (4) ensuring the existence of a periodic state in a system with a NPE 
closely resembles the usual form of such condition. In spite of the apparent analogy, however, there is a concealed 
distinction as regards the actual phenomena which occur in a NPAS as compared with an ordinary nonlinear one. 


In a NPAS, in contrast to ordinary nonlinear systems, the set of the admissible (dimensionless) frequencies 
@ = 2n/N of the required periodic motion is a countable set. So,one knows precisely all the possible frequencies 
of such motions in advance. Therefore, the determination of the frequency of the periodic motion with the aid of 
the Eq. (4) is reduced to finding out whether the occurrence of a given specific system at any of the known countable 
number of (dimensionless) frequencies is possible or not. 


This particular feature is reflected in the formula for the equivalent complex-valued gain of the NPE, in other 
words, in the expression for the fundamental harmonic of the response of an NPE to the signal (2). 


4. Fundamental Harmonic of Response of NPE To Sinusoidal Signal 





A nonlinear pulse element, unlike the ordinary nonlinear element, performs not only the nonlinear conversion 
but also the “pulse* conversion of the input signal (2). Each of these conversions brings about concentration of high- 
frequency components in the spectrum and the combining of both transformations leads to a twofold modification of 
the original spectrum. 


As a result of these transformations, harmonic components of the composite frequencies appear at the output 
of the NPE. The composite components whose frequencies coincide with those of the input signal (2) are superim- 
posed, and in the aggregate they represent the fundamental harmonic of the output signal. 


The “pulse® transformation of the spectrum brings about a repeated translation of the initial spectrum to the 
high-frequency region. Hence the fundamental harmonic at the output of the NPE is the result of superimposing the 
composite components generated by these translations. 


The computations of the fundamental harmonic at the output of a NPE could actually be based on different 
physical features of the investigated pulse elements, and the fact that the composition of harmonics occurs need not 
be taken explicitly into account. Primarily, however, the phenomena are of the same nature. 


We shall now concentrate on a general method of determining the fundamental harmonic at the output of a 
NPE; the above-described pattern of building up the fundamental harmonic will also be illustrated by an example 
with an amplitude NPE. 


The process at the output of a NPE is the outcome of the modulation of the carrier train of instantaneous 
pulses (4). Therefore, the formula for the fundamental harmonic at the output of a NPE can be written directly as 
the result of the modulation of the fundamental harmonic of this train. 


We represent the carrier train of the instantaneous impulses with period T, in the form 


N-—1 


57, (t) = yr yr, (t —s). 


s=0 


The expression under the summation sign is 
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dyr, (t —s) = nF; \! + 22 cos l= G—s)], 
and the complex-valued amplitude of its fundamental harmonic can be written as 
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ied 


2 
a (8) = 77R-e © 


The complex-valued amplitude of the fundamental harmonic at the output of the generalized NPE (with the 
modulation of amplitude, width, position and/or shape of the pulse) can be written in the form 


N—} 


Z, = *- A,. (5) 
s=0 


In the above expression 
2 
A,= a.KpO,, (7 +) a; (s + B,), 


where a, Bs and y, are the dimensionless amplitudes (the sign being taken into account), the displacement and 
the width of the s-th pulse respectively (Fig. 3); Kp > 0, a constant, is its nominal amplitude; e, (j@) is the 
complex-valued spectrum of the s- th impulse. . 


The expression (5) for a NPE which generates rectangular pulses takes the form 


N—1 27 xy 
_ » 2 Ty iy (848, +2 (6) 
Z1= Kp 2 a sinF- e ( bs 


The complex-valued amplitude of the fundamental harmonic at the output of a nonlinear pulse element can 
also be given in the form of a series such that the frequency interpretation of the fundamental harmonic at the out- 
put of an NPE becomes evident. This series can, in the general case, be directly obtained for the amplitude- pulse 
element. 


_ Amplitude NPE We represent the nonlinear amplitude- pulse element (Fig. 5) 
RT OTE OF OR Se AY , MRS ata YG = as a series combination of an ordinary nonlinear element NE, an 
u(t) | v(t) w(t) ) 2) ideal pulse element IPE and a shaping device SD with the frequency 
characteristic Kye. (j). The ideal pulse element generates the 
L--~—~--~~-~~-~~----- 4 instantaneous pulse train such that the areas of the pulses are made 
to be proportional to the values of the input signal v(t) = v(n). 


The response of the nonlinear element to the signal (2) shall 
be written in the form 
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and the train of the instantaneous pulses at the output from the ideal pulse element as 


i 
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Hence, the complex-valued amplitude of the fundamental harmonic is 
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rl co 
W=F, > Bitzi, (7) 
l=—oo 








where ] — afi +N—2E (F )] (A = 0, +1, +2,...), and E(x) denotes the integral part of a number x. 


Direct from (7), the expression for the complex-valued amplitude of the fundamental harmonic at the output of the 
NPE is obtained: 


co 
= > By, (8) 
l=—oo 


1 » 2% \ yp’ 
Bysin = F KpOy (; Bin. 


For the amplitude-nonlinear PE generating rectangular pulses, the expression (8) becomes: 


ay pi ; 
Z,= Kp sin “ > Bin. 


The fact that the expression (5) is identical with (8) in the case of an amplitude-nonlinear PE can be easily 
verified by substituting the actual values of the parameters, a, =v{[s], Bs=U, Ys= 7 in (5). 


This gives 


N-—1 


= 2 As = 7K p? (74 Pe wire 


and the summing with respect to s yields 
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us/2 One can obtain either directly or by substituting the re- 
spective values of as, 8, and ys in (5), the series expressions 
for the first harmonics at the output of a width- or a time- 
| NPE, In contrast to the series (8) whose terms can usually be 
10, 1/2 e expressed by elementary functions, the series for the width- 
or time- modulation represent the fundamental harmonic in 
| | | terms of the Besse] functions of the first kind. 
HITS ST } 7 we The expressions (5) and (8) enable one to draw some 
— general conclusions as regards the fundamental harmonic at 
t-0 vill Ling mn the output of a NPE. 
a oe \ ae It can be seen from (8) that the complex-valued ampli- 
tude of the fundamental harmonic at the output of the ampli- 
vad ———— ane w tude- NPE is up to a constant factor a sum of complex-valued 
1) tye, amplitudes of fundamental harmonics (for a given element) 
and of an infinite number of high-frequency components of the 
7 WaT. Hh} response to the signal (2) of the ordinary nonlinear element 
Su, ' (Fig. 6). The series for the fundamental harmonic at the out- 
12 | H | put of a width or time-NPE are essentially of the same nature 
—+} PSUR ST ON as regards the forming of the fundamental harmonic in the 
; 2a, = 10, latter cases. 
: : : | | | H | | x There is a well-known exception to the above rule, 
— taaiies namely the amplitude- NPE with smooth, polynomial nonlinearities. 
In this case the fundamental harmonic consists only of a finite 
Fig. 6. number of composite components since now there is only a finite 
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number of harmonics in the response of the ordinary nonlinear element. If the characteristic of the nonlinearity is 
written as 
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- 
” - >) ayurt 
=0 


the expression (8) for the complex-valued amplitude of the fundamental harmonic at the output of the nonlinear 
pulse element can be reduced by means of elementary rearrangements to the form 


’ 1 a 
2, = 7, K,01(i 7) 
(ir) nee oi +1 


. pa > 224 (, }1+-1N | -1) ayC2tt+eit+INo 
t— —————__ , 
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The distinctive feature of the nonlinear transformation representing a NPE liesalsoin the response of a NPE 
to a sinusoidal signal; it differs from a response to an ordinary nonlinear element as the former depends not only on 
the amplitude but also on the frequency and the phase of the input signal; moreover, it also depends on the shape of 
the pulse, the repetition rate and the mark-space ratio of the pulse element. 


5. Solution in Frequency Space 





The equivalent complex-valued gain of the NPE can be obtained from the formulas (5), (8) and (2). It follows 
from (5) and (2) that 


N-1 


(CN, ¥, Fo, Dy Ti) = ZS) Ae. (9) 


s=0 
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Fig. 7. 
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For every value of & = 21/N, to the right-hand side of Eq. (4) there corresponds a domain in the amplitude- 
phase plane. The graphicalsolution of Eq. (4) is reduced to finding the points of intersection of the frequency 
characteristic depending on three coordinates, the modulus, the phase and the frequency, with the domains corre- 
sponding to different values of the frequency = 2n/N(N = 2, 3,.. .) (Fig. 7). In actual computations the re- 
spective “projections® are studied on the coordinate planes; this enables one to use the usual amplitude- phase, 
amplitude- frequency and phase- frequency characteristics. 


The determination of periodic states in systems with amplitude or time modulation is conveniently carried 
out by relating the shaping device of the pulse element to the linear part of the system. The calculated frequency 
characteristics now take the form 


K (jo) = Kp®, (7@) K (jo). 
The Eqs. (5), (8) and (9) can now be simplified by the factor 


K,®, (j 3). 


At the output of a pulse element with a relay modulation characteristic the periodic process [that is,the re- 
sponse to signal (2)) can assume only a finite and specified number of different forms for any given frequency 
@ = 21/N whatever the type of modulation. Correspondingly, for each frequency @= 2n/N the amplitude and 
the phase of the fundamental harmonic at the output of the NPE can only assume certain specific values. 


This facilitates the evaluation of an equivalent complex-valued gain of such “relay” pulse elements; this 
amounts primarily to finding the boundaries of domains of values of the gain corresponding to different values of the 
fundamental harmonic at the output foreach @ = 2n/N in the plane of the amplitude-phase characteristic. By 
carrying out the corresponding transformations, the expressions for the fundamental harmonic can be brought to such 
a form that the equations of these boundaries can be determined directly. 


6. Use of the Describing-Function Method in the Evaluation of Periodic States 


in NPAS. 





In applying the describing-function method in the analysis of nonlinear automatic systems one presupposes that 
the fundamental harmonic of the periodic process at the output of the linear part LP considerably exceeds the higher- 
frequency components of the process. For this assumption to be valid, it is necessary that the LP behaves as a low- 
pass filter [13, 14]. 


In the majority of cases of practical importance the linear part of the nonlinear pulse system satisfies this 
condition and the describing-function approach can be used as effectively as in the case of ordinary nonlinear 
systems (see, for example, the comparison of the numerical and experimental data in (8, 10, 11)). 


The describing functions can be more successfully applied in NPAS by using shaping devices in pulse systems. 
These, as it is known, behave as low-pass filters and are especially used to smooth out pulse control signals. 


With the aid of the shaping signals the behavior of the “pulse* and the “continuous” part of an NPAS can in 
principle be always coordinated in such a way that the level of the high-frequency components at the output of the 
LP of a nonlinear pulse system will only differ a little from the one at the output of an ordinary nonlinear system. 
Other things being equal, the higher the repetition rate, the mark-space ratio of the pulse element as well as the 
order of the shaping device,the lower is the level of high-frequency components at the output of the LP of a non- 
linear pulse system. 


Restricting the passband of the LP of pulse-reproducing systems (in compliance with the requirements of the 
Kotel'nikov theorem on the sampled-data representation of a continuous signal [15]) provides favorable ground for 
the describing function being efficiently applied in the analysis of pulse servo systems, discrete-to-continuous data 
converters, devices with extrapolators, programmed control systems, etc. 
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SUMMARY 


The main features relating to the use of the describing function in NPAS have been discussed already; it was 
shown that this method remained also valid for systems with various types of pulse modulation as well as with ar- 
bitrary nonlinearities in their pulse elements. 





The obtained results can serve as a basis in actual computations to determine approximately the periodic 
states in various cases, When evaluating them, one only meets with technical difficulties of a purely computational 
character. These difficulties can in many cases be overcome relatively easily by an efficient use of the specific 
properties of the pulse elements under investigation. Here an NPAS with pulse elements which have relay modula- 
tion characteristics could serve as an example. 


ee nee oe 


The formulas for the fundamental harmonic were obtained with only slight restrictions imposed, and they are 
therefore applicable to practically any nonlinear pulse elements. However, one should not rule out the possibility 
that just because of their universality, these expressions can prove in certain specific cases to be less useful than 
others obtained easily by making use of the special characteristics of the pulse element under consideration, From 
the experience gained in the use of the describing function in NPAS (see Introduction) one knows that only rarely 
does one come across these “easy” ways. They should not nonetheless, be ignored when tackling specific problems. 





The linear part of the system was described by its frequency characteristic K (j@) as is usual in the describing- 
function method. It was shown recently by Ya. Z. Tsypkin [16] that by using in Eq. (4) the pulse frequency character- 
istic K *(j@) instead of K (jd) one is able in some cases to obtain a better approximation, and in very simple cases 
(N = 2, 4) even an exact solution. 
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LINEAR CONTROL OF LINEARLY-ASYMMETRICAL OBJECTS 
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pp. 1493-1497, November, 1961 

Original article submitted January 28, 1961 


The paper studies relay control of an object which is characterized by different pure lag 
times or different transfer functions for the closed-loop and open-loop states of the relay element. 
Such objects are called linearly- asymmetrical by the author. The problem of determining the 
oscillation period of the system is investigated. 


A number of heat power engineering control objects have the feature that the heating and cooling time con- 
stants for the object are not identical. In a number of cases the pure lag times for heating and cooling the object 
ate also different. Not to take into account this asymmetry in the computations (as S. F. Chistyakov pointed out 
[2]) means to obtain computed data which is radically different from the experimental data. Objects with the 
properties indicated above shall be called linearly- asymmetrical (LAO). We shall study systems containing an 
asymmetrical relay element (a two- position controller) and one of two types of LAO: 1) an object in which only the 
pure lag values differ; 2) an object with no pure lag but with different transfer functions for the on and off states 
of the relay. Since in such systems the relay element (the relay section of the system) usually has comparatively 
large moving masses (a contactor, a slide- valve with a solenoid or a motor drive, etc.) it is necessary to verify the 
system over the minimal period or for the maximum oscillation frequency of the system so that the indicated 
frequency does not exceed the allowable frequency for the relay section of the system. 


1. A Linearly-Asymmetrical Object with Different Values of Pure Lag 





A furnace whose heat transfer is regulated by steam-valve control of the entire unit or of its individual sections 
can serve as an example of such an object. Here the pure lag time when the steam valve is closed will exceed the 
pure lag time which occurs when the steam valve is open; the difference will be equal to a value that is approxi- 
mately equivalent to the condensation time for the steam which is trapped in the furnace. 


We shall denote the transfer function of the object when lag is neglected by W(p) and the pure lag time for 
closure and opening of the relay shall be denoted by ry and rp. 


When these systems are studied, we shall make use of the theory developed by Ya. Z. Tsypkin. To the extent 
possible we shall choose the same notation as that used in [1). 


Figure 1a illustrates an unsymmetrical characteristic of a relay element. The notation here is the following: 
Yo is the value at the input of the part of the object that is not switched off. For the case of controlling an incom- 
plete flow y, # 0(3], + xo is the control range, 2K, is the magnitude of that part of the object which is switched 
off. Figure 1b shows the variation of the control signal y (1) as a function of time, Here T is the period of the 
oscillations, y, is the relative time for which the relay element (RE) is switched on, z(t) is the controlled quantity, 
f, is the controller setting. 


The quantity [z(t) - f,] is applied to the input of the RE. It is evident that z(t) arises as a response of the 
LAO to an infinite number of pulses that are of height 2K, with a period T, a relative “on” duration ,, and a con- 
stant effect y,». Here during the period from 0 to t, the transfer function of the system contains a constant pure lag 
T and during the period from t, to y it contains a value of pure lag equal to ry. 


We shall assume that the system has identical pure lag times equal to rp. 
Then 
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z(t) = 2K, {h(t —%) + >) kh (KT’—1,) —h(K ~1’) tT) o<t<rT) 
K=1 


and 
z(t) 2h, |- h(t—t,) + h(7'7T’ —%) 
— Sh(KT" —1,) —h(K 4-7) T — 
K=! 
(T’ <tr’, 
where Z(t) = z(t) —yW (0). 
Originating from these values, we determine the characteristic of the relay system 
’ a" & ” -~ re 
J(T") = — an z,(7")— jz,(T") 
and 


, 


- A A Pe , oe 
It (PT) = — x24 (iT) —F2, (117). 
Taking into account the fact that z(T) =z (0), h(0) = 0, for 


t<%, —h(t—t,) = h(0), 





n oo 
' al 
h(t.) — h (tm) || > cy, (e¥'e — ef ‘m) and >) a® = - ~ for a< 1 
v=l K=1 
(the values c,,, Py and n are the same as in (1)), we obtain 
f es ePullT ee) Pyle To) 
z(7 ) = 2K, {3} Cy, a ae | 
v=] | —€ 


and 





~ = mae PYUT'=%) PUY Te) 
(7',7") =2K,[3 Cy, Cara ee Pe. Lit snk. ° ||. 


From this, we have 


dhe : Py (LT’—*.) Py (lv) 7" —to] 
("= 2K |e + ¢, $= 
v=] — eé 





and 


7s 





; ~ 4 PyUT =) Py LW IT sh 
~ , On ae ho) e — 
Serre an, (—e 4 Spe [rorr—w_ er arin 7 


v= l—e 


(the values c’, cy, are the same as in (1)). 
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The characteristic of the relay system is 








P - P, (IT’—*,) Py (U—*’)T’--=] 
Tr. ( \) a” P ois ° 
ie Ana , . ’ > 
ee eee ST a. PLT — 
vol | — ¢€ 
Pere PUT’ —. * (l—v')T’—= 
F on \! ‘ e vf Bi i ae mt ' )T o] } 
5 Ce RT Sey > are | 
v=} i—e 
and 
n P, (IT’—=,) P, (I+) T’ — <0] 
T’ ‘ , : : : i ee’ , Pig i o 
J , (T’) =_— — A Be Cc + bs Cc ef vl o) ee, . a 
Tt c} — Ye Pi 
v=1 | P 
— PUT’ —=,) P. ((l+y")T’ —- 
+ ° 2x y' P, 609 ; Te) 1 e 2 ~€ vl | 
PT a fw |! We be Sue. Gpe eon > one 
v=] | -€ 


We shall plot J,(T") and i, s(T*) in the complex plane for a series of values yj while varying T* from 0 to ow. 
On the same plane we shall draw horizontal straight lines which intersect the imaginary axis at the points 
inn’ f vt Xo + %W (0)and — f y— %o + YoW (0). The points at which these straight lines intersect the root loci 
J, (T*) and J, »(T*) make it possible to plot two curves in the coordinates y* and T'. The point at which the latter 


two curves intersect yields the actual values of y* and T*. 


y Analogous computations should be performed while assuming that both pure 
lag times are equal to 7,; in this way we obtain y * and T*. 





The actual system will behave as a symmetrical system in which the pure 
| c lag is equal to tr, when the relay element is “on*; i.e., under these conditions 
+ the behavior of the actual system will correspond to the behavior of a symmetrical 
40 system with a pure lag equal to r, over a period of time y‘T*'. In Fig. 1b this 
time is the interval from 0 to tp. 
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le When the relay element is switched off the behavior of the actual system 
will correspond to the behavior of a symmetrical system with a pure lag equal to 
Ty over a period (1 - y")T® (i.e., over a period during which the relay element 
of the symmetrical system is switched off). In Fig. 1b this time will be the in- 
terval from tz to ty. 
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It is evident that the total period T = y ‘T* + (1 - y")T" and the relative 





= 


time for which the real object is switched on is given by y = a : 





< 


Determining a series of values for y*, T* and y", T® for various values of 
f,, and plotting the graph [y'T* + (1 - y")T"] =f(f,), it is possible to determine 
the minimal period T and the controller setting for which the period will be 
Fig. 1. minimal. For the furnace this corresponds to determining the minimal period and 
the furnace-intake air temperature for which the period will be minimal, From the 
same graph it is possible to determine the period for any value of the controller setting, 





_ 

fF. o. — — — + 
ctr - — - + 
SP h. ne ee oe 
rr --= 
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2. A Linearly- Asymmetrical Object with No Lag Transfer Functions 





An example of such a system can be seen in a system with an object for which the relay element connects 
and disconnects a linear correcting element simultaneously with the basic circuit. 


The block diagram for such a system is shown in Fig. 2. In the figure LO is a linear object, CS is the correct- 
ing section; RE is the relay element. The remaining notation is the same as above. From the block diagram it is 
evident that in the case where the contacts of RE are closed the transfer function for the linear section of the system 
is determined by the transfer functions of the linear object and the correcting section; in the case where the contacts 
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of RE are open the transfer function of the linear section of the object is equal to the transfer function of the linear 
object (this follows since the output of the correcting section can be disconnected from the input of the linear object). 


In that case we should perform the same operations as in 
Y% section 1 in our computations; we assume / =1 and +r =0, In 
2ft) order to obtain y * and T* it should be assumed that the transfer 
function of the LAO is constant and equal to the transfer function 
when the relay element is closed. In order to determine y* 
and T* the transfer function is made equal to the transfer 






























































= i ol function of the LAO for an open relay element. 
en's The value of the relay system characteristic (for deter- 
cs mining y* and T*) is 
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Fig. 2. omnt i—e” 
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The subscript of P,, denotes that these quantities apply to the transfer function of the object under conditions 
for which the relay element is closed, 


From the analogous characteristic (but for the other transfer function) we determine y “ and T*. 


It is not difficult to plot the transient response if we know the poles of the transfer function and take into 
account the fact that over the time intervals 0 to ty, ty to tg, tz to ts, and ts to t, the transfer functions and the 
values of pure lag are different. 


A system with an object similar to that described in section 1 was computed, The computed data and the 
results obtained from simulating the system coincided. These results also coincided quite closely with the results 
of tests made on an actual object. 


APPENDIX 
The object which we experimented on was a mine furnace with steam valve control. 


A theoretical derivation and an experimental check were made of 
the proposition that the unit is an aperiodic section of the first order with 
pure lag, Here the pure lag r» when the valve is opened is determined 
solely by the time required for the air to averse the distance from the 
furnace to the spot where a low-inertia resistance thermometer is mounted, 
The pure lag +, which arises when the valve blocks the steam pipe exceeds 
T» by the amount of time required for the valve to block the steam pipe 
and the time required for condensation of the steam trapped in the furnace. 


In the case under study the object has the following specifications: 
The time constant of the object is T, = 200 sec; r,=18 sec; ry = 28 sec; 
Xo = 0.045; @ ma, is the maximum increment of the output parameter, 
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Control is achieved on the basis of full intake (i.e., yp = 0). The system was theoretically designed and simulated 
on an *MN-7* simulator with a relay attachment. The variation of the external temperature (or the variation of 
the average control level) corresponds in our case to a variation in fy. 


The dependence of the oscillatory period of the system on the variation of fy is given in Fig. 3 (f, is given 
in fractions of the maximum increment of the output parameter). The solid line is used to denote the computed 
curve.* The crosses denote the results of simulation; the circles denote the values obtained as a result of the ex- 
periment. | 


A value for the period that is somewhat too high was obtained during the experiment; from our point of view 
this is due to the fact that the operation and release of the steam pipe valve mechanism did not occur at the instant 
that the controller contacts touched because of the oxidation of those contacts. An additional time was required to 
produce a certain pressure on the contacts. 




















- | as | 13 | 12 | 2/3 | 56 
a | 0.17 lo 3 | 0.47 | 0.61 | 0.776 
Pa. bs | 0.18 | 0 33 | 0.47 | 0.62 | 0.78 
T ( 

exp | 9-187 | 0.43 | 0.525 | 0.585 | 





The table cites data for the comparative time in accordance with the computation, simulation and experiment, 
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ON PROGRAMMING PROBLEMS FOR DIGITAL DIFFERENTIAL ANALYZERS 
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The paper cites a method of scaling problems for solution on digital differential analyzers 
(DDA) that operate with a binary code and use the binary system for representing increments. 
Such a method is used at present in solving radio-engineering and electrical communication 
problems on the universal digital computer “Integral” in the M, A. Bonch-Bruevich Leningrad 
Electrical Engineering Communications Institute. 


Programming problems for digital differential analyzers reduces to formulating the functional network for 
connecting the integrators and-to computing the scales of the variables. These two stages of programming reduce 
the problem to a form which makes it possible to program it and introduce it into the machine, The method of 
formulating functional networks and programming various mathematical functions on a DDA almost completely 
coincides with the method of formulating functional networks that has been used for analog computers. We should, 
however, note that in solving many problems on DDA wide use is made of special integrator modes which do not 
apply in the case of analog devices. The specific features of a DDA are also manifested in scaling the variables. 
The purpose of scaling is the representation of the variables of all the integrators in a definite scale which permits 
the DDA to solve a problem with a specified degree of accuracy and at a definite rate in accordance with its 
capabilities. The specific result of scaling must be the determination of the length® of all the integrators and a 
notation for the initial conditions which corresponds to their capabilities. 


It is expedient to begin the computation of the scales by choosing the step Ax and the scale S, for the in- 
dependent variable; these quantities are related by the expression 





and determine, as we know, the accuracy and the speed of solution. It is evident that the integration error € over 
certain variation interval for the independent variable does not exceed the quantity Ax | f* (x)| 3, where | f °"G)| nay 
is the maximum value of the first derivative of the function in the specified interval. Therefore the quantities Ax 

and S, can in a first approximation be determined from the specified value of ¢€ on the basis of the relationships 


Ar = epiibesedanle . (1) 
iF) ax 
S, = log, Kf Siler (2) 


Note that in addition to the integration error the digital integrator introduces a whole series of other errors 
associated basically with the binary system which is used to represent the increments (1-3). Therefore, it is ex- 
pedient to reduce the value of the step Ax obtained according to formula (1) by several binary orders of magnitude. 


After choosing the scale of the independent variable the scales of the output variables for all the integrators can 
be computed from the formula [1] 


* The length of the integrator is assumed to mean the number of computation places in the register Y of the in- 
tegrator which are used in solving the specified problem. 
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S, = S,—m, (3) 


where m is the number of binary places assigned for the storage of the integer part of the integrand function y in 
the specified integrator. The quantity m for each integrator is chosen for the condition 


2" > 4! max (4) 


Substituting the value of m for each integrator into formula (3) and thus traversing the entire contour of the 
functional network, it is possible to compute the scales of all integrand functions Sy (and then the lengths N of the 
integrators) from the formula 


N=m-+S,. (5) 


It would be possible to terminate the computation of the scales with this since the values found for N and m 
fully determine the writing of the initial conditions in all the integrators. However, certain special features of the 
functional networks for connecting integrators in the DDA impose additional requirements on the computation of 
the scales and give rise to the necessity of correcting the latter. The correction of the scales is required, in parti- 
cular, for summation of several variables whose scales have been found to be different on the basis of a preliminary 
computation. Assume, for example, that a function f , (t) is formed at the output of the integrator 1 (Fig. 1); 
assume further that the maximum rate of change |f j(t) |; ,ax of this function corresponds to the maximum (i.e., 

a value equal to unity) rate of overfilling the specified integrator. In that case the register Y of the integrator will 
be completely filled (m, = m, ip). * Such a utilization of the integrator register is most rational, since a) the 
available range over which the rate of filling can vary (from zero to unity [3]) is fully utilized, and b) the maxi- 
mum possible accuracy of the solution is achieved for the specified integration step, In fact, since m = Mpjin under 
these conditions, it follows that the quantity S; =S, - m_ proves to be the maximum possible value for specified S,. 
However, such a utilization of the register is not always possible. Assume that the function f,(t) which must be 
combined with f, (t) is formed at the output of the integrator 2; here | f} (t)| max >? lf2 (0 max’ in the addition 
it is necessary to use the condition Sz. =S,, ; therefore the maximum rate of change f,(t) cannot correspond to the 
maximum (i.e., a value equal to unity) rate of filling the integrator 2, since this rate corresponds to | f, (t)| max ° 
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scale dy m =! f.Scale Sy, 
min |’ 
m-m, f(t) 4.Scale Sy 
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— F(t *: 5 in the 
—t— -—_ scale 
: _ oe 
min |2 
as ee 
, scale 5,, §. Scale 5,25; 
a Scale 53*% 
Fig. 1. Correction of scale for summation of Fig. 2. Block diagram for raising the scales. 
variables, 


In order to achieve Sz, it is necessary to place the integrator 2 under conditions where its register Y is not 
completely filled (i. “ it fe not be used rationally). Such a mode of operation is reached for my, > Mgmin- In 
other words, in the integrator 2 it is necessary to lower the scale of the output variable; this is achieved by increasing 
mg i.e., by increasing the number of places that are reserved for the integer portion of f3(t). Lowering the scale is 
the simplest case of correcting scales and is used both for summing certain variables and for reducing the overflow 


* Min ‘8 the minimum integer for which the equation holds. 
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rate; for example, it is used when a tracking integrator is overloaded. The addition is achieved by increasing the 
scale of the variable. 


Figure 2 cites an example of such an increase in scale by means of an additional (scale) integrator. Assume 
that it is necessary to sum the variables f,, f, and fs on the scale S,, and assume S;,= S;,-= S,, and 
S;,--Sx—1, since my min =1. Increasing the scale S f, by one can be achieved by means of an additional in- 
tegrator scale 2 analogous to theintegrator 1. Summing the overflows of both integrators in the scale S$, we obtain 
S + =S,; here division of the variables f,; by two does not occur since two increments in the scale S, are equivalent 
to one increment in the scale S, - 1. 


The special feature of functional networks for connecting integrators in DDA, as was noted above, is the ex- 
tensive utilization of special integrator modes which differ radically from the conventional integration mode. Most 
frequently a tracking integrator is used which performs the functions of an adder which adds several variables (1, 4). 
This summing is realized not by the addition of the increments of the variables but by the addition of the rates of 
these variables for constant scales. If we use pip y’ Pings - - -» Pinn © denote the rates of the variables at the input 
of a tracking integrator, and p,,,, to denote the overflow rate for the tracking integrator proper, then the following 
relationships will be valid: 


Pn (Pin, + Ping + +--+ + Pin,) for | Ping! 4 Lies Ping | <1 (6) 


and 


Pout = +1 for Pin. Ping 00 Pin, 4. (7) 


It is evident that in performing the inequality (6) «..e sum of the increments at the output of the tracking in- 
tegrator per unit time will be equal to the sum of the increments of the variables at its input. In the contrary case, 
the tracking integrator will cease to respond to any changes in the sum of the variables applied to its input until the 
register Y of the integrator returns to the original state with the maximum positive or maximum negative number. 
Such a tracking integrator mode (an overload mode) is completely impermissible. Note that after overload occurs 
the output rate remains equal to + 1 even in the case where inequality (6) is satisfied. Therefore, it is evident that 
correct programming provides for satisfaction of inequality (6) at the input of the tracking integrator over the entire 
period of solution of the problem. In order to estimate the maximum value of the sum appearing in the left side of 
inequality (6) we perform a rough computation that reduces to finding the upper limit for the rates of the variable at 


the input to the tracking integrator. This computation consists of the following. For any integrator, is is possible to 
write 


) <ip r 
out ' max Pin'max max (8) 


if 
where rp is the coefficient for maximum filling of the register Y in the integrator: 


RAP 
Tmax = es (9) 
Based on the maximum rate of the independent variable of the problem (if this is the machine rate it is equal 
to unity) it is possible to use formulas (8) and (9) to compute the values of r yay and |P for each integrator; 


out | a 
then the possibility of overloading the tracking integrator can be estimated from inequality @. 


One of the methods for eliminating overloads of the tracking integrator consists of the above-mentioned 
lowering of the scale for the overflow inputs; this is analogous to reducing the rate of the variables. 


Another method is based on using special network diagrams for connecting the integrators [4]. An example 
of the practical realization of one of these is shown in Fig. 3. Here Au, Av, Aw are quantum increments of certain 
variables which are summed by the tracking integrator. The idea governing the operation of the networks resides 
in the fact that each increment at the output has a weight which is twice as great as its value at the input; because 
of this the network can produce the sum without overloading. Note that in contrast to the conventional tracking in- 
tegrator the specified network reduces the scale of the increments by one. We should also mention that the network 
in Fig. 3 (and similar networks) are used expediently only in DDA of the parallel type, since only in this case ts 
complete symmetry in the operation of both integrators achieved, 
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Thus, in accordance with the analysis above it is possible to recommend the following sequence for programming 
problems on the DDA. 


1, From the specified differential equation we formulate a preliminary functional network for connecting the 
integrators. 


2. By analyzing the differential equation on the basis of inequality (4) we determine the values of m for all 
the integrators. If it is not possible or practical to perform such an analysis, then the values of m are specified 
arbitrarily. 


3. From the specified accuracy of the solution we choose the scale of the independent variable in the problem 
while taking into account the specific features of the functional network. The expediency of the chosen scale Sx 
can be estimated from the point of view of the time required to solve the problem; this time is computed from the 
formula t = A2°x /k, where k is the speed of operation for the machine (the number of iterations per second), and 
4 is the integration interval. 


4. From formula (3) we compute the scales S,, S, and Sy for all the integrators. We perform the necessary 
correction of scales (i.e,, the lowering of scales by means of a corresponding increase in the value of m in the in- 
tegrators, and an increase of the scales on the basis of special diagrams for connecting the integrators). 


5. In the case where tracking integrators are used the relationships (9), (8) and (6) are used to estimate the 
possibility that such integrators will be overloaded; then measures are adopted to eliminate overloading. The 
corresponding modifications are made in the diagram. 


6. The lengths of all the integrators are computed from formula (5). 


1. A operative diagram for connecting the integrators is formulated. All of the computed quantities are 
indicated in the diagram (Fig. 4). 


s 



























































Trackin 
inte grat 
fe 
4z fr: - 
—| Trackin . 
{Eee 2 au sdorber 
Fig. 3. Programming network which avoids Fig. 4. Notation used for the integrator in the 
overloading the tracking integrator. operative diagram (n is the number of the 


integrator). 


8. The problem is programmed on the machine in accordance with the operative diagram. If the values of 
m for all the integrators are chosen correctly, it follows that the problem will be solved completely without the 
Sietiniaaii of the programmer. However, if the values of m are chosen arbitrarily, then certain integrators may 
overflow during the process of solution and the machine must be stopped to change the scales. 


In conclusion, we shall study examples which illustrate the proposed programming method. 
Example 1. We shall formulate the program for solving a nonlinear equation that describes the operation of 
a rectifier with a filter consisting of one reactive element: 


dy : 

dz + ¥ =/ (sinz —y): 
j/(sinz—y)=sinz—y for sinz—y>0, 
/(sinz—y)=0 for sinz—y<0 
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for the initial conditions 


y(0)=0, yy’ (0) == 0. 


It is known that | < 1.Assume further that the absolute error of computation must not exceed € = 0,001, 
'dz | max 

1. The functional network diagram for connecting the integrators is given in Fig. 5. In integrator 1 we com- 
pute the quantity y= /fy'dx; integrators 2, 3 perform the function sin x; integrators 4, 5 are tracking integrators and 
form the sum sin x +(- y); the integrator 6 is the null-organ which transmits control by limiting the function 
sin x - y; integrators 7, 8 serve to divide the function 2 (sin x - y) by two. Division by two by increasing the scale 


of the increment cannot be employed, since the increments arriving at the input of the integrator 1 must be in an 
identical scale. 


2. From the maximum values of the integrands we 
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Fig. 5. Operative network diagram for Fig. 6. Operative network diagram for 
programming the nonlinear differential programming the formulation of the function 
equation. u =xe~ %. 


3. The scale of the independent variable of the problem is 


1 
ct, 06 Dey 
S,> loge 0.001 


and we choose S, = 15. 
4. Further, we determine 


S,, = 5, —m = 15-1 = 14, 
So, =S,— im = 15—1= 14, 
S.,4 = S,—m = 15—1= 14, 
Syg = Sig = Sy'= 14, Sy = Sy, 
— He Oe 


== Sig = Uh, Sia = Sey = 1h, Sue 


5. For the tracking integrators 4,5 we verify inequality (6): 


i 


1 
“ r = odes ca =e 4, 
Pin 1 max pe in 3 max’ 3max xt! 2 : 
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Thus condition (6) is satisfied. 
6. We determine the lengths of the integrators: 
Ni =m +S, =1+ 14 = 15, 


Na = N= m+ Sy = 1+ 14 = 15, 
Ne=m+S,,=0+1=1. 


Example 2. We shall formulate the program for computing the function u = xe %* for a = 0.3 in the interval 
0=<x<=10. 


1, The functional network for connecting the integrators is shown in Fig. 6. Integrators 1, 2 form the function 
v =e” &%; integrators 3, 4 form the product xv; the unknown function u accumulates in the integrator 5. 


The initial and maximum values of the functions in the integrands are equal to 


yi (0) = v’ (0) = — 0.3, y (0) = — 0.3 = const, ys (0) == 0, ye (0) = v (0) = 1, 
Ys (0) — 0, 
1 Imax = %3+ 1Yslmax = 1% | %lmax = ! 
2. Based on the maximum value of the functions in the integrands we choose m, = — 1*, mg= — 1, ms =4, 


mM, = 1, 
3. We choose S, = 14. 


4. We determine 
Sa = Sq, — ™ = 14+ 1 = 15, 
Si. = Sig — ™ = S,, — ms = 15 
+ 1= 16, S,, = S,, — m=S,, 
—mg=15—4=11, Sig = Sug 
—m, = 14 — 1 = 13. 


It is evident from the computation that to satisfy the equation Sz, = Sz, in the integrator 4 it is necessary to 
lower the overflow scale by two units. Therefore we choose my = 3. 


Now 
Sy = 14—3=11=S,,, Si, = Sz = 16, Syg = S, = 14, Sy = Sy = 15, 
Sug = Szq = Sy = Al 
6. Therefore, 
Ni =m + Sy =—1+ 16 = 15, Ns = mg + Syg = 4+ 14 = 18, 


Ne = m+ Sy =3 +4 15 = 18. 
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ad 


* m =— 1 denotes, in which the maximum value of the function in the integrand does not exceed 27'. 
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ON THE PROBABILITY CHARACTERISTICS OF COMPONENT RELIABILITY 


S. E. Rostkovskaya 

(Moscow) 

Translated from Avtomatika i Telemekhanika, Vol. 22, No. 11, 
pp. 1504-1512, November, 1961 ' 

Original article submitted February 22, 1961 


A method for analyzing the probability characteristics of reliability, establishing their 
relationship to the law of distribution of the component's parameters, is described. Examples 
of the analysis of such characteristics obtained experimentally are given. 


In the commonly-accepted quantitative definition of reliability the criterion is the probability of perfect 
operation of an apparatus (or its components) during a prescribed period of time. It has been pointed out repeatedly 
in [1-3] that it is expedient to classify faults with respect to their different effects and the measures which can be 
used to reduce their number. Among the most frequent causes of element failures we mention the following: 


1) mechanical breakdown; 
2) electrical breakdown; 
3) irreversible physical changes resulting in changes in the component's parameters. 


We shall assume that the operating conditions are such that direct mechanical damage to the components 
does not occur. Damage is only observed as a result of electrical breakdown. With a sufficient margin of electrical 
stability the number of failures in the initial period of time is negligible and their time distribution is subject to the 
law of infrequent phenomena (the Poisson distribution). In the course of time, as a result of chemical changes in the 


material and surface phenomena, changes in the electrical parameters of the components take place and the break- 
down voltage decreases, 


Let us consider the relationship between the law of distribution of failures, at the time when the results of the 


above-mentioned processes begin to make their appearance, and the laws of distribution of the element parameters, 
including the breakdown voltage, as functions of time. 


We shall show how to determine the law of distribution in the time interval T— T,, during which the value 


of the random function p (representing a certain parameter of the component) reaches a given limiting value pjjm.- 
Below we shall term the time interval T — T, the service life. 


We first consider the method for determining the time interval T— T, at the end of which the value of p 
exceeds the value pjj,, for 50% of a single given component type, taken in sufficient number. For this, we assume 
that the law of variation of the median p meq = f (t) is known, and that the function is monotonic along the entire 
axis and continuous at all points (Fig. 1). We note that this characteristic is also the characteristic of Tpeq 

= f (pip) with substitution of the coordinates, pj;,,, in place of p and 

4 Tmed in place of t. The length of time T = Tmmeq is determined from the 
ZA point of intersection of the curve p meq = f () with the straight line 
CD plotted parallel to the t-axis at a distance p jmp, comesponding to the 
of Pred permissible range of variation of the parameter p. 


7 i LL Then for the triangle with sides Ap and At =AT, for which two 
LG Xp Fun 
a 





vertices are located on the curve pmeqg= f (t), with sides parallel to 
the coordinate axes, we can write [4] 














» (Tmea) AT = @ (Peg) Ap, 


Fig. 1. 
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where ¢ (p is the probability density corresponding to the value pmed, ¥ (T meq) is the probability density 
corresponding to the service life T,,.4. In the limit, as A p ~ 0, we obtain 


qPmea 


» (7 mea) = P (Pirea) —— 








Similarly, we can plot the curve of values p; = f (t), satisfying the condition 


Pi 


F (p,) = | (9) dp = const 


—oo 


and at a distance At from p meq = f(t) at the ordinate pj = p)jjm- 


For the curve p; = f(t) corresponding to some constant value of the distribution function, i.e., variation of a 
quantile in time, we obtain 


dp; 
~ (Ti) = O(P,) ae - (1) 


Thus, the probability density of the service life of the component T = Tj is equal to the probability density of 
the parameter value pj = pjj,, at this instant of time multiplied by the derivative of the curve corresponding to the 
constant quantile pj = f(t) at this point. 


As an example, let us consider variation of the median of the random quantity p with time according to the 
law 


t 
t 


Pmea =i—e 


We put the change of p; corresponding to a certain other quantile, taking place according to the law 


t 


A=P,+1i—e *, 
where Py is a realization of the random quantity, whose distribution law we assume normal and invariant with time 


i ini all 
=e 20* 
 (P,,) Vi’ 


To obtain the law of distribution of service life according to (1), it is necessary to determine the derivative 








om 1 - + 
_—. 2 
We obtain the required law in the form 
' ey ‘ Tj ( 2, 2) 
Pe eee oe. A> agile 


16 V dn : 


where Py and Tj correspond to the points of intersection of the curve pj= f (t) with the straight line p = pjim. 
i.e. sal Gi 
° Py = Pum —1 +6 ?. 
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Let us construct graphically the service life distribution curve. For this we plot on the graph (Fig. 2) the 
distribution curve of the parameter (a normal distribution is assumed), the characteristic of time variation of the 
median and the curve of variation of pj; (t), corresponding to the condition of constant probability 


Pmea = / (t) and Pea + Py = J; (t)- 


The point of intersection of the straight line, plotted at a distance p);,,, from the axis of abscissae, with the 
curve representing the variation of a certain quantile as a function of time, determines the service life, whose 


d 
probability density can be defined as p(7';) = @ (p;) | pa . In practice, according to available data, the distribu- 


tion law of the parameters of some component gauged by individual operations is close to normal. On the contrary, 
in sorting completed parts it approaches equiprobable. 


Aside from this, the parameter scatter varies with 

time and its distribution also varies, and therefore it is 
te necessary to investigate the variation of the distribution 
law with time. However, under identical external con- 
ditions and load, there is no basis for supposing that a 
substantial deviation of the distribution law from normal 
occurs, It is more likely to assume that the parameter 
scatter increases while the law remains constant. 


Y , 1 T The distribution of time intervals up to the instants 
(?) Uy WL : at which the parameter passes beyond the operating 
¥ YY Pim tolerance depends on the magnitude of the latter. An 
Z 
Y 


P=pyt!-& 


¢l 1) Ped =1-@ “ 











| increase of the permitted magnitude of parameter deviation 
corresponds to an increase of median service life.* The 
service lives of elements for which the realizations of the 
random values of the parameters are close to the limiting 
values admitted in the circuit may be substantially iess 











Fig. 2. than Tneq (T min in Fig. 3). The scatter of service lives 
depends on the shape of the function pj, = f (t) satisfying 
the conditions of constant probability level. 

R : 
rae” 
-tke 
Ye) (4) Pmea*'-? = 





\ 
\ 
\\ 
\\ 
\\ 
\ 
\ 
| 
| 
| 
| 
\ 
| 
| 
! 
! 
| 
! 
am 


_— 
n= 
—. 














a 








| 
2 


Fig. 3. 


* The normal distribution law is understood. 
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Let us consider a variation of the quantity Py with time, i.e., we assume that the curves p{(t) are not parallel 
to each other. We assume in addition that with time the dispersion of the random quantity corresponding to this 
parameter increases. Let the variation of the quantile as a function of time take place according to the law 


where Prao is the initial value, The change of dispersion of the random quantity in this case is defined by the 
expression 

o? = (1 + kt)* 
(see Appendix). 


We now determine the distribution of service lives corresponding to a normal distribution of the parameter p, 
under condition (2) 


(ey = Pmed)” (1+-4t)?? 2 
1 - 207 (1+41)? _— i Py [20° 


V 2s (1 + kt V 2x0 (1 + kt) 








~ (p,) = 





Substituting the latter expression in (1) we obtain the distribution law for service life 


% |S 
i — <3 1 ar. 
7h on ail 20 
¥(7) 314 kT) Vin 





These examples illustrate the method of determining the differential law of component service life distribu- 
tion. 


To determine the relative number of elements whose parameters exceed the permitted limits at a given in- 
stant, we employ the integral function 


Tj 
Q= \ w(ryar (3) 
or, on the basis of (1) 
Q@=1— | o(P)dp, (3°) 


where the limit of integration p; is defined by the condition p ;(t) = p jim. 
Employing this expression, we can also define the reliability characteristic 


at 


P=1—Q= \ o(p)ap. 


—oo 


The laws of distribution and time variation of the parameters also present interest since the exponential 
reliability law 








p, 


3) 


3") 





Employing (3) and (1), we obtain 


a | 


Q@=1—\ ear 
T; $y * r 
- \ H(T)dT=1t—e *, 


whence 


- 


T 
‘mo od {-— 
(7) =e * or PP) = —e ’ 





Fig. 4. Consequently, an exponential law of service life distribution 


y (T) can be obtained, for example, with an equiprobable distribution of the element parameters in the initial (or 
some other) instant of time and an exponential law of variation of the quantile values of the element parameters 
with time (Fig. 4). 


T brs. 


, 





Fig. 5. 


Let us use the above method to analyze the curve of the time distribution of failures for a certain type of 
ceramic capacitor (Fig. 5) [5]. In the graph,the time scale along the axis of ordinates is logarithmic, and the scale 
along the axis of abscissae for the number of failures corresponds to the probability of their appearance with a 
normal distribution law Q (log T). 


In [5] an empirical relationship is also given for the service life of the capacitors as a function of applied 
voltage with Q=const., i.e., for a certain number of permitted failures: 


T.=T, | =) oF (4) 


where n= 3, z= rere is a function of ambient temperature, T, and T, are the service lives of the 


Capacitors at temperatures 6, and 0x, respectively. 
























The curves of distribution of the number of failures were obtained at a constant temperature 150°C, hence 
Z =comnst, 


1% On the basis of the data available, using expression (1), we can 

ws find analytically or graphically the distribution law of capacitor- 
breakdown potential at a certain instant of time, for example, at t 
=T,. For this, we define from (4) the function 


— z 
afte nV ua of. (5) 


According to the graph of Fig. 5, in place of y (T) we shall 
consider y (log T). From (5) we find 






Q2 


05 
1 
2 


5 


logu, =logu, + (logl, —log7;) 








+ + zlog2, (6) 
whence / (108%) — _.1 and according to (1) 

d (log T..) 3 

1 \| 4 (loge) 

lo => l i» ee eee . 
Fig. 6. Pq (0B Mx) = Yo(log?'s) d (logT,) 
1 

Consequently, Pq (log Ux) = > Yq llog7’). 


Varying T, according to different values of Q,we find the law of distribution of the logarithms of capacitor- 
breakdown potential, which is also normal. 


Q% 


r 


400 500 600 700800 1000 1200 





Fig. 7. 
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In order to plot the distribution ¢ (log u) from Fig. 5, in which we represent the integral curve of the total 


number of failures with increase of experiment duration, we find the differential curve and correspondingly (log T) neg 
(Fig. 6). Similarly to (6), with @ = const, we obtain 


log ux —logu, = + [(og7") mea — (logl").}. 


Plotting this line and the straight line 


log u, —logu, = + [log7’)mea — (log7’)» + @), 


where 9 is the quantity by which the straight line is shifted along the log T axis with change of Q, for a given instant 
of time, for example T = 1000 hrs, we plot the distribution curve of log u,. 


With a more careful consideration of the graph of Fig. 5 we note that the straight lines Q= f (log T) plotted 
on it are not parallel, and therefore in the formula for the dependence of service life on potential for various values 
of the number of failures the exponent n will have different values. Aside from this, it is unlikely that there will 
be greater scattering of log T at a voltage of 900 V than at 600 V, although the duration of time to failure in the 
latter case is greater than in the former (Fig. 7). It should be remarked that the points obtained when testing the 
capacitors at 900 V have a greater scatter with respect to the straight line than at 600 V. It appears more correct 
to plot the initial portion of the line for 900 V parallel to the line for 600 V up to 60% failures, where a time in- 
terval appears during which no failures are observed, The remaining section of the graph should remain without 
change. Carrying out this construction, we find that the distribution of log T and log u in this case do not corre- 


spond as a whole to the normal law, but may be approximated on individual sections by curves of this form (differ- 
ing in mean values and dispersions (Fig. 8). 
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SUMMARY 


The law of distribution of the number of component failures (service lives), defined by the time of deviation 


of a parameter beyond the limits permitted in a circuit, is related to the law of distribution of the parameter in the 
following way. 
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The probability density of the number of failures is defined by the product of the probability density of the 
parameter and the derivative of a certain auxiliary function at a point whose ordinate at the given instant is equal 
to the limiting permissible value of the parameter. The family of auxiliary functions represents the dependence of 
the parameter on time, for each preassigned constant value of the quantile. Determination of the mutual dependence 
of the distribution laws permits approximate prediction of the probability characteristic of the reliability on the basis 
of the parameter distribution law at the initial time and the forms of time variation of the parameter median or 


some other quantile, 


To increase the reliability of circuit operation, it is necessary to select those components for which the 
production tolerances on the parameters, taking into account reversible and irreversible changes, do not exceed the 


permitted variation of the parameter in the circuit. 
The author expresses appreciation to B. I. Filippovichfor valuable advice in discussing the article. 


APPENDIX 
We define the dispersion as 


DP yg (1 + kt)= Mp?_, (1 + kt)? — [Mpy_, (1 -+ ke)}?. 


Since the factor (1 + kt) is constant at a given instant of time for the distribution of the random quantity, it 
may be placed outside the sign of mathematical expectation 


Dey (1 -i- kt) = (1 + kt)? [Mp?_, — (Mp,_,)*I. 
The expression in square brackets is the dispersion of the random quantity p at the initial time. Consequently, 


DPi—9 (1 + kt) as (1 -+ kt)? Dp; 2 (1 t kt)? o?. 


LITERATURE CITED 


1. N.G,. Bruevich, On the reliability and precision of automatic production [in Russian] Izv. AN SSSR, Otd. 
tekhn. i, Energetika i avtomatika No. 4 (1959). 

2. G.V.Druzhinin, On failures in automatic systems with aging or wear of elements [in Russian] Izv. AN SSSR, 
Otd. tekhn, i., Energetika i avtomatika No. 4 (1960). 

3. S.-M. Kuznetsov, Probability of damage to automatic system elements [in Russian] Avtomatika i telemekhanika 
19, 11 (1958). 

4. I. V. Dunin-Barkovskly and N. V. Smirnov, Theory of probability and mathematical statistics in engineering 
{in Russian) (Gostekhizdat, 1953). 

5. Hi-Reliability Ceramicons, IRE Transactions on Reliability and Quality Control, PGRQC-11 (August, 1957). 





All abbreviations of periodicals in the above bibliography are letter-by-letter transliter- 
ations of the abbreviations as given in the original Russian journal. Some or all of this peri- 
odical literature may well be available in English translation. A complete list of the cover-to- 
cover English translations appears at the back of this issue. 


1372 





of 
ence 


he 


tly, 


R, 


anika 





ON CERTAIN PROPERTIES OF SECOND-HARMONIC MAGNETIC 
MODULATORS WITH SINGLE-PHASE AND TWO- PHASE 
SUPPLY CIRCUITS 


M. A. Rakov and L, A. Sinitskii 

(Lvov) 

Translated from Avtomatika i Telemekhanika, Vol. 22, No. 11, 
pp. 1513-1520, November, 1961 

Original article submitted January 21, 1961 


The paper analyzes two types of second-harmonic modulators: modulators with single- 
phase and two-phase supply circuits. Based on a comparison of the gains and time constants 
for both types of modulators, the advantage of the two-phase modulator over the single-phase 
one is demonstrated. 


Second harmonic magnetic modulators can have single-phase or two-phase supply of the excitation windings 
(1). The extensive use of modulators with two-phase excitation was faciliated by the appearance of sources which 
produced two voltages in quadrature [2, 3]. 


For two-phase supply of the modulators the control circuit does not contain a second-harmonic current. This 
makes it possible to obtain a high gain without using a filter in the control circuit, and to achieve a lower time 
constant, 


However, we should not mechanically identify a two-phase modulator with a single-phase modulator in a 
mode where the even current harmonics in the control circuits are suppressed. The shorting of harmonics of order 
4n in the control circuit has a substantial effect on the operation of a two-phase modulator, Therefore, it is of 
interest to attempt a study and comparison of the properties of two-phase modulators. 


In a two-phase modulator there is no necessity of employing an additional filter in the excitation circuit when 
the modulator is fed from a quadrature phase shift network [3]; thus the modulator can continue to operate when 
appreciable fluctuations of the supply frequency occur, Therefore, it is natural to compare such a modulator with 
a single-phase modulator that has a choke in the signal circuit, since only such a “simplest” filter can be assumed 
to be sufficiently independent of frequency. 


The effect of the inductive filter on the gain and time constant of a single-phase modulator has been estimated 
previously for assumptions which introduced appreciable errors in the computations, Therefore, for an objective 
comparison of two types of modulators it proved necessary for us to investigate a single-phase modulator with an 
inductive filter. This investigation is of independent interest, since it provides recommendations for the selection 
of such a filter and permits a more accurate determination of the time constant. 


The analysis which took into account the higher harmonics permitted both an objective comparison of the 
modulators and a quantitative estimate of the advantages which were achieved when a two-phase modulator was 
fed with square voltages as proposed in [2]. Both the analysis and comparison were made for no-load. 


A Modulator with Single-Phase Supply 


The circuit of a modulator with single-phase supply is shown in Fig. 1. We shall write the differential 
equation of the signal circuit as 





; : d 
LS + riz + Sw, = (By — Bu) = Us. (1) 
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Here S$ is the cross-sectional area of the modulator cores. The flux densities in cores I and Il are functions 
f (1) of the sum and difference of the field intensities H, and H, that are produced by the excitation and control 
currents. We shall assume that the control signal (and therefore 
i=] sinwt the field intensity H,) are so small that the Taylor expansion of 
the differences in fluxes Bj — Bj; can be written as 











By — By = 2f'(H,) Hy. (2) 





The dependence of the dynamic permeability on the field 
intensity is approximated by an expression corresponding to the 
arc tangent approximation of the function B = f (H): 








’ 1 
By=/ (4) eh, See 


Here ji, is the initial magnetic permeability, and H, is the value of intensity corresponding to pg = 3 /2 Ho 
For a sinusoidal excitation current 


| 
Mg Ho u \’ , (4) 





1 + (a7 sin et 


8 ; 


where | is the average length of a line of force. After introducing the dimensionless parameters 


wl, Wala ese r ° wel Ly 
Reig ts: maa tat te oe EE ART: 





2 
where L, = 2p, 2 is the initial inductance of the signal winding, Eq. (1) will become 


aa J, Poe. mee ee 
We + the + 8 how erg rh, (5) 


where (2) and (4) are taken into account. 


In determining the gain K,, of the modulator the voltage drop across the resistance r can be assumed equal 
to zero for all harmonics 1,(h,) with the exception of the dc component h{?), such an assumption is possible in 
view of the large inductive reactance of the filter choke. Equation (5) under these conditions can be separated 
into two equations: 


d ; hy on 
ag (te +4; ry =; ) 0, (6) 
no ods ho. (7) 


Integrating (6) and solving it relative to ha, we obtain 


b= C[1—B;— oe], 


where C is the integration constant, 


hi 
} 2 
B == 2 ’ a ; 2 . (9) 
h2 h? 
1+A+—- +A + a 








ae 








From (8) we find the second harmonic of the dimensionless field intensity, After determining the constant C 
from (7) and (8), we obtain 


DP sii oe AE 
Pm te (10) 


Since there are no emf harmonics in the signal circuit, the variable components of the fluxes in the choke 
and in the control winding are equal. Therefore, the second-harmonic flux associated with the control winding is 
equal to 


y? — LIA? = Linas Vira 2 == 2B hp, 008 2, (a) 


H,! 
where /], = ' 
Ws 





The emf in the second harmonic of the output winding is equal to 


a  “edp®  SolwU; 14~-yi-a .. . 
ee a rae @—Vinai sin 20. (12) 








The modulator gain is Ky, = =e? ) t/U2, where E(2) denotes the amplitude value of the second-harmonic emf; 
on the basis of (12) this coefficient rig oak to 


r i—Vi—a? 
ui = 8v i (8—Vi—a) ; (13) 





oLlwe 


Wer 





where v = 


As the filter inductance increases, the gain increases; for L + o it reaches a maximum value corresponding 
to absence of even current harmonics in the signal circuit: 


h2 
6 deat 
~s So (14) 


Ku 1mex= 8v| === — 1 
ee Le a 


The graphs in Fig. 2 determine the dependence of the gain on the value of the filter inductance L. 


We shall now estimate the time constant of a modulator with an inductive filter. Since the transient response 
is caused by the establishment of a dc current component, we should study the full Eq. (5) in this case, For con- 


venience in analysis we shall introduce the new variable /, — he é 
1 +4, sin*.§ 





Then (5) becomes 


(1424 & cos 20) q+ ly (148m) 4 ht, (sin 20— 1 cos 20) | b— thy. (15) 


Equation (15) is a linear differential equation, and therefore the nature of the transient response is completely 
determined by the form of the solution for the corresponding homogeneous equation: 
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i ¥ +a \ 
, 0+ arc tan(|/ i tane) 
b= b, — 


(16) 





(1 —acos 269)% +1 


n= . ° 
2(1-+ A) + Az, 





where 


The angle 6 enters into the exponent of Eq. (16) in nonlinear fashion; therefore, rigorously speaking, we 
cannot assume that the process involved in reaching a steady state is exponential. However, if the duration of the 
transient response is greater than the period of the external emf (i.e., the interval of variation for 9 is much greater 
than 1 ) it is permissible to write the approximate re- 
lationship arc tan(k tan 9) = 9, since this equation is satis- 
fied exactly for @ = nw /2, where n is an integer. Then 
the time constant will be equal to 
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Fig. 2, 
Cerf -42— Ditibiedletmalienilia 7: 
: ( V (2+ m2, 42-2) — nt 


The expression for the time constant is conveniently represented in the form 


L+Leq 


r 


Here Leg = pla is the equivalent inductance of the control windings: 


(17) 











2 
V (240,42 FY m2 


Relationship (17) is shown graphically in Fig. 3. The decrease in the coefficient p as hy increases is an 
easily explainable physical fact. An interesting and not so obvious fact is the decrease in p as the choke inductance 
L decreases, This fact must be taken into account in determining the modulator time constants for small values of 


filter inductance. 
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A Modulator With Two-Phase Supply 





The circuit for a modulator supplied from a two-phase source is shown in Fig. 4. The equation for voltage 
equilibrium in the signal circuit is 


dB dB dB dB 
wyS ( in 5% t+ — = uta ri, = U3. (18) 











The voltage at the output of the modulator is 


Tour = —w2S a (Bim — Bum— Bin + Buy). (19) 


Assume that the excitation currents for the halves M and N of the modulator are respectively equal to ij 
=I sin wt and iy =I), cos wt. 


t=], sinwt 























Having made use of Eqs. (2) and (3), we obtain 











AH 
Bim — Bum = 2po = 
1+ (FPsiner) 
Hy, 
Bin — Buw = 2po aI ; (20) 
1+ WT coon | 


Having substituted (20) into (18) and (19) and made use of the dimensionless parameters introduced in the 
preceding section, we obtain 





2(2+h2,) rc *li. |e eae (21) 
Yo zy. — j 093*20 


where Yo= 





oly 


Oour == — 25H, UpSwh 





at cos + (22) 


m rn) r —— 
ae =) _ = cos? 26 


In order to simplify the solution of (21) we replace the expression in brackets by a new variable y. We obtain 
the equation 


mY + (1 + d008.48) y =a, (23) 
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2+ him hin . ho 
’ = ,a= J 
(hi, + 8hF, + 8) Yo ime 48 °° AS + Bh? 48 
Equation (22) is written as 





where m=8 








d 
U at an 2poH we Sahin Fr (y cos 26). (24) 


out 


The linear differential first-order Eq. (23) can be reduced to quadratures; however, under these conditions we 
obtain integrals which cannot be expressed in terms of elementary functions. For this reason, it is expedient to use 
approximate methods to solve this equation. We shall seek a solution for (23) in the form 


y = Ay + A, cos 46 + B, sin 40. (25) 


Having substituted (25) into (23) and set the coefficients of the trigonometric sums on the left and right sides 
of the equation equal to one another, we find 


2 
Ne = its z= 4, A, = 6 a, B, = mb a. 


tm? — > [a 1+m—* 











Based on the derived expression for y we find the second-harmonic voltage at the modulator output 


= Bh;,V Us 
out AS + 8248 44 mi _ © 


2 





U [(2 + 2m? +- 6) sin 20 — md cos 26], 





The gain of a modulator with a two-phase supply is written as 











4 2 3 2 
Ku =—— ors torte (26) 
hm + +8 1 +m? — — 
m 
on the basis of this equation. 
If we. > 1, which is usually the case, it is possible to neglect the terms containing 6 in (26): 
8 
Ku. = . ° (27) 
i+ +8 
m 


The time constant for a modulator with two-phase supply can be found by solving the homogeneous equation 
corresponding to Eq. (23). The solution of this equation is written as 


~+ (¢— {sin o) 


Y = Yoe 


Since 6 = wt, the time constant for the circuit is equal to + = 


e |5 
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des 


26) 


(27) 


ion 





Expressing m in terms of the modulator parameters, we find 


24-h? 
g Le | m 


so , (28) 
r hf + 8h2,+8 


t= 





From this the equivalent modulator inductance is 


2+ h2, 
Leg = ~9 - * men 4 a 
8 + 8h? + hf 





Comparison Between Modulators With Single-Phase and Two-Phase Supply Circuits 





In order to judge the expediency of one of the two possible methods of supplying magnetic second-harmonic 
modulators, it is necessary to compare the gains and time constants of both networks. Comparing (13) and (27), it 
is not difficult to prove that the dependence of the gains on the structural parameters (this relationship is determined 
by the value of the parameter v) is identical for both types of modulators, The values of v for identical windings 
and core data coincide in both cases if we neglect the choke resistance. 


The difference between the gains is caused solely by terms which are characterized by the mode of operation 
(the parameter h,,). The effect of the mode of operation for a single-phase modulator with an ideal filter and for 
a two-phase modulator is determined in accordance with (13) and (27) by the multipliers 


hrm 
F, (hm) eae (lim) 
1m) = | 57> F 72 2 Vm pe nn 
Vi 1. f2 R h2 1 8 
tm hos 3 & 


The graphs of the functions F,; and F, are shown in Fig, 5. From the graphs it is evident that the maximum 
gain for two-phase supply is shifted toward smaller values of hy, than those corresponding to the maximum gain for 
single-phase supply. Moreover, the rate at which the gain decreases as h,, increases is greater in a two-phase 
modulator than in a single-phase modulator. The sharp drop in the gain of a two-phase modulator as the excitation 
increases is evidently due to the growth in harmonics that have a number which is the multiple of four and have 

their paths closed through the signal source, 
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However, in a single-phase modulator, it is necessary to introduce a filter that appreciably worsens the 
modulator characteristics in order to achieve a mode close to the ideal mode; the characteristics of the modulator 
deteriorate both due to the voltage drop across the filter and due to an increase in the time constant, 
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As the criterion for an objective estimate of modulator quality we may take the ratio between the time con- 
stants of the two networks when their gains are equal. The relationships derived above permit such an estimate to 
be made easily. For this purpose we find the ratio ¢ = F,/F, for different values of h,, from the graph in Fig. 5; 
this ratio is equal to the ratio between the values of gain for the two-phase and single-phase modulators, Based on 
the value of ¢,we used Fig. 2 to find the value of the filter inductance L for which the single-phase modulator will 
have the same gain as the two-phase modulator. Finally, based on the value found for L, it is not difficult to deter- 
mine the value of the time constant for the single-phase modulator and the ratio o = r,/r__ for the time constants 
of the networks under study. 


The graph for the function o = f(h,,) is shown in Fig. 6. From this graph it follows that for all modes of 
operation the modulator with two-phase supply offers substantial advantages over modulators with single-phase supply; 
this follows since for identical values of gain the two-phase modulator has appreciably smaller time constants. 


The advantages of the two-phase modulator become even greater when its excitation windings are fed from 
square-wave voltage sources (as proposed in [2]), since in such a case there are no even harmonics above the second 
in the signal and output windings, 
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DYNAMIC CHARACTERISTICS OF ELECTROMAGNETIC 
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Based on an analysis of the electromagnetic powder clutch (EPC) as an element having 
distributed parameters,the dynamics are examined, and its time and frequency characteristics 
as well as transfer functions are determined, 


An engineering design method, universal dynamic characteristics, and a sample cal- 
culation of the EPC dynamics are presented, 


Electromagnetic powder clutches (EPC) are finding ever greater use in automatic speed control systems, in 
servomechanisms with a large number of reverses, in fast-acting servo systems, etc, [1-3]. In connection with these 
it has become necessary to study EPC characteristics in detail. 


EPC static characteristics have been studied by O. N. Tatur, As for dynamic characteristics, they have been 
studied only experimentally in isolated, special cases [3] in order to get essential quantitative data for the computa- 
tion of specific automatic-control systems. Moreover, the EPC time characteristics were assumed to be exponential. 
Such a description of dynamic characteristics for an EPC, which is an element having distributed parameters, is 
correct only for certain conditions (when there is appreciable lag in the excitation circuit) that usually impair its 
dynamic behavior. 


In the present work results are reported on theoretical and experimental studies of the transient responses in 
EPCs from which their dynamic characteristics have been determined, The investigation was carried out under the 
following universal assumptions: 


1) The mechanical characteristic of the clutch M (n) is absolutely fixed, i,e., 
M (n) = const; (1) 
2) The moment of the clutch is proportional to the magnetic flux 
M=k,®. (2) 
Note that conditions (1) and (2) are not compulsory for the solution of the present problem, but only make the 


analytical study somewhat easier, Subsequently, these results are corrected for the case when the moment is depicted 
by the relation 


M = ky®", (3) 


where 1 <n <1.5, depending on the construction of the clutch, 
1. General Expression for the Dynamic Characteristics 
Definitions. When evaluating the dynamic properties of an EPC we will consider the following characteristics: 


1) The natural dynamic characteristic of the EPC, which is the reaction of the clutch to a unit voltage step 
U =1(t): 





Dy = Dx(t), (4) 
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My= Mx(t), (4a) 


where ©), and My are the magnetic flux and the rotational moment of the EPC (in absolute or relative units); 
2) The driven dynamic characteristic* of an EPC, which is the reaction of the clutch to a unit current step 
i=1(0: 
@® = M(t), (5) 
M = M (t). (5a) 


The need to introduce the magnetic-flux characteristics, which do not appear as an output quantity for the 
EPC, is explained by the substantial effect of the magnetic circuit on the dynamics of the EPC, 
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Actually, on the strength of Eq. (3) the EPC dynamic properties are determined by the electromagnetic processes, 
i.e., expressions (4) and (5). Expression (5) is the principal one since fast-acting clutches are used mainly in the 
current supply condition, The characteristic of (4) is subsidiary, but it turns out to be useful for selecting an EPC 
control circuit. 


Peculiarities of the Electromagnetic Processes in an EPC. Operational Form of Representation for the Dynamic 
Characteristics, When an EPC is connected to the control voltage (Fig. 1), a transient response occurs in it which 
is determined by the lag due to the excitation circuit and the massive magnetic circuit. 


The effect of the excitation circuit lag for a reasonable choice of control circuit parameters is secondary 
while, on the other hand, the lag of the massive magnetic circuit plays an essential role. 


A general solution to the problem of finding the electromagnetic transient responses in a clutch with a massive 
magnetic circuit is given in [4]. 


In order to find the dynamic characteristics of an EPC on the basis of this solution, it is necessary to develop 
expressions for the operational reluctances f;(p) which are applicable to the magnetic circuit of a powder clutch 
(Fig. 2), Compared with the circuit investigated in [4]** this circuit has the following features, 


1, The presence in the working gap 6 = const, of a ferromagnetic layer having a relative permeability yp, 
in the order of 4 to 8, 


2, A considerable thickness of the inside pole, which is represented by the solid cylinder: 
t=e (0<e<1) 


re 


* In the future whenever the subscript N is absent the driven characteristic is implied. 
** See expression (22) in [4]. 
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The first feature makes it necessary to allow for y; when figuring the equivalent permeability [4] 


eS R+PL 7: 
< Ri+ Ry 





Ue = (6) 


The second feature requires that the operational reluctance f,(p) of a thick-walled cylinder be determined 
(see Appendix I), The appropriate expression has the form: 


a) for large values of p 





x (i -- 1+ — 
hie(p) = 7 ; 5) ( toy pr, (1) 
b) for small values of p 
S 1 -- pT 
lec (p) ~ 1 +- (1—8) pT ° (8) 


Calculation of the parameters for the disk portions is carried out in accord with formulas (6) from [4]. 


With the calculation of relations (6), (7) and (8) thus obtained, a general solution (22) from [4], which deter- 
mines the transient response for a step in the voltage on a clutch, gives the natural dynamic characteristic of an 
EPC in operational form, The corresponding time function is represented in the general case by an exponential 
series, 


When a clutch is fed directly from a source having a small internal impedance, 


k k k 
pT >>) %/i(p), >} aifi(p) = >) ui = 14, 
i=1 i=1 i=1 
and (22) from [4] degenerates to the ordinary exponent which establishes the law for the change in current and flux, 
When a clutch is fed through an extremely large resistance Rp » r, or when the internal impedance of the 
source is high, a current supply condition exists.* In this case 


k 
b> afi (p) > pT, 


i=1 





Lo 
where 7), = R + Re and a, are the relative reluctances of the massive parts (fixed parts). 


Expressions (21) and (22) from [4] acquire the form 





i(p) = J, = const, (9) 
@(p) { 
* ala rege’ : (10) 
SY a4; (P) 


i=] 
Expression (10) represents the driven dynamic characteristic of an EPC with respect to current, 


The latter case corresponds to a realistic relationship for the majority of automatic-control systems having an 
EPC, in particular,when an EPC is connected as a load in the output stage of an electronic circuit. 


* A similar condition can also be obtained by connecting into the control circuit driving elements, for example R 


orc, 
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From (22) of [4] and (10) there results the following link between the transfer functions corresponding to the 
natural (4) and the driven (5) dynamic characteristics: 


1 | 
=e ( = ———-... 
© 1 /On- PTo ° Dy | / ® + pl ’ 


2. Calculated Transfer Function and Time Characteristic of an EPC 





Current Supply Condition, Using the general expression for the dynamic characteristic of an EPC (10) and 
taking into account that for large values of p in accordance with [4] and (7) 


h(p) =uV p, (11) 


we get the calculated transfer function of an EPC in the high-frequency region for the current supply condition 
(T, & 0): 


@M(p) 1 





=—. (12) 
 tVP 
The summation of 1; gives.the following result: 
. ite) -- 
t= > 4 (1 + é,) — VT;. (13) 





i=1 


For disk parts it must be assumed that &; =0 and e€,; =1, since these coefficients have significance only for 
cylinders, 


To find the calculated transfer function in the low-frequency region there can be utilized a representation of 
(10) in the form of a rational fraction or a linear averaging [4]. 


However, a quadratic averaging gives results which are more favorable for plotting the characteristics, 


Let us return to the parameter 7 of formula (12) and assume that the time constants of the massive parts are 
identical, i.e., Tj =T,. Then,squaring (13), we obtain 


4 si 
Te = 12 k (14) 


{= 








Thus quadratic averaging is equivalent to the replacement of a compound magnetic circuit in the clutch by 
some equivalent massive element having a parameter r identical with the real magnetic circuit, i.e., an identical 
characteristic in the high-frequency region. 


In the case of a single massive part (k = 1) quadratic and linear averaging give the same results, as would 
follow from a comparison of (39) out of [4] and (14) for k =1, Tj = Tg, a4 =1. 


On the basis of (14) the calculated transfer function of an EPC can be written for low frequencies, i.e., for 


small values of p, by using the first term of the series for 1/ f,.(p) as the equivalent of the magnetic circuit having 
a time constant T,: 





oe) Bne B "7 
tess ts Y = T+ pT, + (1 — Be), 


n=1 


om (15) 
Do 1 -- pT, ‘ 
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where 


k 


B= >} Bi. (16) 
i=1 
For a uniform cylinder 6 = sar = 0.69. For a flat disk 
g=5%—0.81 


In real clutch designs 6, = 0.77 to 0.8, in which the larger values of 6 ¢ are related to clutches with a larger 
diameter size, 


We will define the time characteristic at high frequencies as the original of the function (12) 


bes (17a) 
Do Vix tT 





Introducing the symbols OM; (t) = ® (t) a= —— ;———~ ie 
0 


we obtain 
®; (i) = =F aVi. (17) 
Va 


The time characteristic at low frequencies is the original of the function (15): 
®; (t) = 1 — Be. (18) 


Solving (17) and (18) together we get the boundary point for the transition from (17) to (18): 


O' (tg) 20.5; tg = 7... (19) 


=~ 0.97, practical calculations can be performed according to the approximate 





a 
(1 + e;) 
Since 5) ai(1 + &:) —5— 
i=1 
formula 


in which the error does not exceed 3% for the possible range of values of €. 


When making accurate calculations to compare different modifications of the magnetic circuit geometry, the 
terms (1+ €;) and (1+ €;) must be computed for given values of €; according to [17], and tg in order to find the 
intersection of (17) and (18), 


Thus, at high frequencies the time characteristic of an EPC (17) is represented by a parabola with an exponent 
of 0,5, and at low frequencies by an exponential function (18) containing a constant component (1 — 8), approxi- 
mately equal to 0.2, 


The parameters of the functions 1, 6¢, and T, which have been derived are completely determined by the 
magnetic circuit geometry and by the ratio p/p of its material that is representative for the current supply condition, 
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Voltage Supply Condition, From (22) of [4] and (15) we have 


DP) _ i 
Do cit tao i op pT. " 
To 
ifd—f)eT, +? 








which gives after transformations 








rf?) _ 1+(1—R) PT 
 ~TFPM +h) +P UBT,” (21) 
Similarly from (5) and (12), we get 
a se (22) 
Dp Tt Vp + pT» 


Expressions (21) and (22) are the transfer functions of an EPC for the voltage supply condition at low and high 
frequencies, respectively. 
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Fig. 3, 


The time characteristic corresponding to the function (21) will have the form 


D, = 1 — Aer! — Bert, 








(23) 
where 
_— —(¢@+iI+VP4+1,29+1 1 
Pie = 0.49 _ (24) 
a= 4 is the electrical lag coefficient of the EPC excitation circuit, and 
e 
(5/T, +7) 1 (5/To+ Ps) 1 ( 
SS _ en SS Se ——— . 25) 
" Pi(Pi— Pr) QT,’ . P2(P2— Pr) IT, 


The original of (22) is expressed by means of the special functions erf x ("probability integral"), which is 
somewhat complicated to evaluate, However, it is not needed here since for a weakly-driven EPC it is sufficient 
to use expression (23), 
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When selecting calculating formulas the q coefficient must be evaluated, For q < 0.1 it is advisable to 
utilize the transfer functions (12) and (15) and the time characteristics (17) and (18), while for q > 0.1 use (21) and 
(23), respectively. 


In particular, for the large values q = 5, expression (23) takes the form 


q t 
NWO _,_ i tbe air. (23a) 
q+1 


-3, = 
Generally, for a 15 to 20-multiple drive, i.e., when 


R, +e +r 


r 


kp = -= 15 to 20 


a current supply condition exists, For special drive cases employing capacitors, nonlinear resistances, etc,, it is 
satisfactory to have kp = 10, 


Disconnection of an EPC by cutting off the Current, In 
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0 rT] relay and pulse automatic-control systems a method of disconnect- 
2 i aa LL 1a( 38) _|$10 wf) ing the clutch often used is to cut off the current of an output 
| mrs: ta eae are irs | || Re | stage. The corresponding time characteristics are easily obtained 
5 | Pet formula (34) \ | A. by subtracting (17) and (18) from unity: 
al co caskcantealh cnn t-) ; ie 
| ‘Dp ier 1 ans +2 } t, 

p> lla sexcedians xe ‘ ine 1 dist” ~ Vx (26) 
Rea WANNA UE WS 2, Oishi ts 
| ~ ee tt i 2disé Ve (27) 
by 1. | N Universal characteristics for (17), (18), (26) and (27) have 
— ae a been plotted in Fig, 3. 
20 

20 log Miu.db In order to find the time scale when calculating a particular 


EPC, one need only compute T, (see Appendix ID), 


Fig. 4. 3. Frequency Characteristics of an EPC 





From what has been presented above it is clear that the most complete advantage of the dynamic properties 
of an EPC is attained under a current supply condition, Therefore, we will consider the frequency characteristics of 
an EPC expressly in this condition. The plot of an amplitude-frequency characteristic in the case of an exponential 
transfer function (21) obviously does not need special consideration, 


Having substituted p = jw in (12) and (15) we get the amplitude- phase characteristics in the form (the sub- 
script on T, is omitted): 
(jo) a 
Do V of 
2 (jo) _ 1+(1—8) jo (29) 
— 1 -+- j@T' . 


(cos + — jsin =. ¥ (28) 





The amplitude and phase frequency characteristics at high frequencies have the form 


W,\@) = ie be She 28 (30) 


0, (@) = -- + = const; (31) 


at low frequencies 
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9. (w) = arctan| -- 


(33) 


On the basis of (30) and (32) we will write down expressions for the logarithmic frequency characteristics 
(LFC) at high and low frequencies, respectively: 


20logW, (w)] = — [4+ 10logwT], (34) 
20log[W, (w)] = 10log 1 + (1 — 8)? w*7'*] — 1010g(1 + w*7?). (35) 


A plot of the curves for (34) and (35) is given in Fig. 4. As seen from the figure, in the range from 0 to 7 db 
the curves differ by not more than 1.2 db (at the point wT = 1), 


We will recall that for high frequencies |wT = 2| the characteristic is found directly from (34). This straight 
line gives the LFC of a clutch over the whole range of frequencies with an error not over 1.2 db, at the point wT =1, 
Thus, for an EPC the concepts of the asymptotic and the actual frequency characteristics are the same and the error 


of 1,2 db at low frequencies is appreciably less than the usual error for an LFC of an inductive section at wT =1 
(3 db). 


The pass-band of an EPC, limited to a level of 10-14 db is 4-10 units of relative frequency i.e., the “frequency 
cutoff” w..= 4/T for 10 db and w.=10/T for 14 db, The LFC of an exponential transfer function gives, as is 
well known, a significantly smaller pass-band. 


Taking the Assumptions into Account, Form of LFC when M=kyo". When n # 1 in formula (3), the time 
characteristics for the moment can be found from the characteristics for the current by raising the latter to the power 
n, When necessary, the transfer functions with respect to moment can be determined by finding representations of 
the functions M (t)/Mg, or they can be obtained directly from the transfer functions $ (p)/@, according to the 


principles in (5). 
In the prevailing case for n = 1 the current and moment characteristics coincide, 


We obtain the frequency characteristic for an EPC when n > 1 by raising the complex quantity (28) to the 
power 0: 





W, (o) = — —, (36) 
(@T)2 


0, (@) =— a = const. (37) 


From (36) we get an expression for the LFC 


20log W, (w) = — [20nloga + 10nlog wT). (38) 


The origin of the LFC is located at the point 


. a] 1 
logWi (o)=0, of =—, 


the coordinates of which are independent of n. Asn increases the LFC rotates clockwise in the direction of larger 
values of 6, changing the slope from 10 to 10 n db/dec, 


For a given value of n the LFC is universal and it is only necessary to compute Te. 


4. Experimental Data 





An experimental determination of the parameters characterizing the dynamic properties of an EPC is possible 
only under those conditions for which the choice of the transfer function is predetermined by experimental conditions. 
Obviously the most suitable cases are for q s 00, and q w 0, 
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In the first case the time constant T, = Ly/r, is found, which may be replaced by measurements of the clutch 
flux and current in the static condition for a known value of ry. 


For the second case, it is convenient to disconnect the clutch with a fast-acting switch (for example, an elec- 
tron tube) thus ensuring the conditions t = 0, 1 = 0 and ZR =, 


In addition, curves of 4j,..(t) and Mgigc (t) are taken on an oscillograph, and Te is found in accord with the 
specification of (17), or + in accord with (17a). An oscillograph record of current is taken in conformance with [4]. 


For measuring the moment the clutch rotor is brought to a stop by means of an inertialess dynamometer of the 
tenslometric type. An experimental curve® of 4j,,.(t) for one of the clutches studied is presented in Fig, 3, Since 
M = © for the particular clutch, the curves of 6 (t) and M (t) coincide. 


We note that the time and frequency characteristics taken experimentally for intermediate values of q, 1.e., 
corresponding to formula (23), can be used only for a graphical estimate of the system. To find from them the 
necessary calculation parameters is difficult. 


SUMMARY 
1, The electromagnetic powder clutch represents a system having distributed parameters, and its dynamic 


characteristics are described in the general case by transcendental transfer functions, 


2. Depending on the class of the drive, the EPC time characteristic is described by an exponential curve of 
the second order (q > 0,1), or by a parabola with an exponent of 0.5 (q < 0,1) in conjunction with an exponential 
curve of the first order which specifies the process for large values of t. 


3, The most advantageous dynamic properties of an EPC are obtained for the current drive condition in which 
the LFC has the smallest slope (about 3 db in an octave), 


4, Owing to the linearity of an LFC the design of an automatic-control system with fast-acting EPCs can be 
satisfactorily carried out by frequency methods, 


5. The relationships obtained permit an engineering calculation of the dynamic EPC characteristics to be 
made with a minimum expenditure of time and an acceptable error (see Appendix I), and also allow systems with 
EPCs to be studied in a general way by analytical methods, 


APPENDIX I 


Operational Reluctance of a Massive Hollow Cylinder 





We will introduce the symbols: J (x), J, (x) are Bessel Functions of the first order, No (x) and N, (x) are Bessel 


Functions of the second order, k = rae is the damping characteristic [6], A, is the selectivity parameter 
47 

("eigenvalve") [6], 6, is the amplitude coefficient of a flux component in the series (6), (1 + €) is an adjusting 

multiplier (see [7]), € = 14/r, is the relative thickness of the cylinder, xp is a root of the transcendental equation 


resulting from the application of the boundary conditions [8], and Cp is a constant of integration, 


For a linear treatment of the problem, the application of Maxwell's equations to the case of a unit step in the 
field (Fig. 5) leads to an expression for the induction B (r,t) which is well known as the thermal transfer equation 


dB eB, 1 aB . 
re (Ta + T Or ) ‘ G-» 
The solution of Eq. (I-1), which is given in [8] has the form in our symbols 


too 


J VA,,"2) 
Bi(r,) = Don exp (— kA,d E (V>”) _ y : 
i 


NV dP) 


* The experimental work was carried out at the driving gear laboratory of the experimental research institute for 
metal-cutting machinery (ENIMS) with the assistance of D, D. Ll"ichev. 


No WV | : (i- 2) 
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Using the method which is set forth in detail in [6] it is possible to calculate the magnetic flux, numerically 
equal to the permeance, by integrating over the cross section. We will give the final result for a positive step in 
the field 


@® 


foe) 
t 
ae - 2 (1— B,e7*n', (I- 3) 


my Ni (Wyn) 
An (ari) (Pa—ri) N2(Vi,r1) — N2 (Vigra) 


—_ z. 4 am 
vi, =— = 
Tey © 


(1-4) 





B, = 








The series in (1-3) is analogous to the series in (7) of [4], In agreement with [4] we obtain the operational 
reluctance for small values of p by neglecting the terms which are quickly damped (n > 1) and writing the re- 


ciprocal of the expression for permeance (I-3): 


























i 
4 4 apie 
f (p) = TPR A 4er (l-5) 
. { ~T-d—Bppr ' : 
5 i+(i— Bs) P= 
pues ‘ : 
where Va; = 2 (1+) 1 _( see {8}), 
} 2 'e--Ti 
H 
| | ten <(1 + 8) < 1.593 |,” 
8 a & i N?(Viar:) 
Nctaaed Th GFA) an) NB (Vann) — 08 (V dara) 
a ra 
Fig. 5. ~ Ia (er) (a—n) (I-6) 
Substituting the value of \, we find 
Ba => SS ae . 
wi+é (i+e) (1-7) 


A curve of (1 + €) as a function of «¢ is given in [7] (pages 269, 325, Fig. 111). 
In the particular case of a uniform cylinder, we have 
mn=0, 1+ = 1533, .,— 2:8 pg. _4 . (i-8) 


2 
"s 


This result was obtained earlier in [6], 


For the calculation of the operational reluctance at large values of p it is necessary to have a solution of 
(I-1) in a closed form, Converting to the operational form of (I-1) and considering p as a parameter we obtain 
Bessel’s equation 


pB (r, p) = k (es ee (1-9) 


B (r, p) = C1 (p) 1o(V-2+) + Cs(p) Nol —% r) (1-10) 
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Making use of the boundary conditions 
B (ra, p) = Bo, D(ri, p) = 0, (I-11) 


we find the constants C,(p) and C,(p), and integrating the induction over the cross section we find the flux 











® (Pp) een 2 
% Ae 
k 
AV =) VE) VF) (VE) 
n(/—2 re] N; (y-2 ri —J; (VY —f ri) No (y-2 re | 
In the particular case of a uniform cylinder 
® (p) 2 (VV -Fn) -13) 





This result was obtained earlier in [6]. 
Like (10) in [4] the last multiplier in (1-12) differs very little from unity for large values of p. 


Omitting the operation for the boundary transformation we will write the operational reluctance for large 
values of p: 


be tau rt 1/ £. (I-14) 


Having multiplied and divided the right-hand part by ¥\, we finally get 


fa (p) = 2 + its V pT. (I-15) 


APPENDIX II 


Calculation of the Time Characteristics for an Electromagnetic Powder Clutch 





Initial Data, Excitation coil: w = 200 turns, r= 0.83 ohm, I, =4 amp, Magnetic circuit: "Yu" steel, p = 
0.35 ohm- mm: m*, Bray = 12,000 gauss, » 4, = 1000. Output stage ("switch"): Type P-208 A transistor, Vinay 
=60 v, Inax=10 amp. Operating gap: 6 =4 mm, y, w 4. Dimensions (see Fig. 2): by =b.y = 2cm, by 
=2.4cm, beg =3.5 cm, fy =3 cm, r2=6.5 cm, fy, = 9 cm, fy, = 11 cm, rj = 4.75 cm, r3=10 cm, lp =5 cm, 


1, Let us determine the reluctances of the parts: 





the operating gap 
8 ° , 
x, = — = 0,75-10 3, up Ry = 0.003; 
bSLPL 

the outer cylinder 

R, = € = 0.04 

a 

1 
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the inner cylinder 


the end disk ("back") 





the disk on the rotor side 


(7, ” ni) 66, 0,047. 


pR,=pR 
” / ci iad 


2. Let us find the relative reluctances 


uw 
R 
4 
I 
te 
= 
> 
i} 
Oo 
to 
_ 


R 
yee = 0.2, @,= no 0.2, a,.= 0.2 
1 > wR, 2 


3. Let us determine the equivalent magnetic permeability 


De +P Lay 


~ 243. 
Dit Ry 





4. We compute the parameter for fast operation and the equivalent constant: 


+ = 112-103) / "e— 0.276, e.= fa — 0.818, e 
ke 


"3 0.462 
p c —. = . 


=" 


1 
id syE [Ody F Od pt boy (1 + Oc) + Abe, (1 + &,)] = 0.43 sec ", 


~ = 0.4231? = 0.077 sec [see (14)and (20)]. 
5. We determine the maximum frequency cut-off at the level of 12 db: 


i 


“a “2 eee 
W; (@) = y= [Wiley] = 86.6 =e 


0.25; 


{ = 
Tt V @ av 
6. We find the time constant of the excitation circuit 


To = 0.4nw* I 


Rn max 


= 0.035 sec 








Here q = T,/Te = 0.45, and therefore to provide the potential high-speed operation during switching it is 
necessary to use a driving network R, and C4. 


The construction of  *(t) [or M* (t)] amounts to the application of a time scale to the axis of abscissas for 
the universal curve of Fig. 3, The scale is applied for the curve of the disconnected clutch. 
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THE POTENTIAL (IDEAL) AND ACTUAL NOISE STABILITY 
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pp. 1533-1542, November, 1961 
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The paper determines the potential (ideal) noise stability of various methods that are used 
for multichannel transmission of telemetering information on the basis of frequency division of 
the channels, A comparison is made between the different methods of transmission in terms of 
potential (ideal) and actual noise stability. A comparison of the various telemetering systems 
with time division of the channels is performed in terms of potential (ideal) and actual noise 
stability in (3). 


INTRODUCTION 





In radiotelemetering systems with frequency division of the channels (RTS-FDC) we must deal with the trans- 
mission of waves rather than with the transmission of individual parameter values (as is the case in radiotelemetering 
systems with time division of the channels), 


V. A.Kotel'nikov estimates the potential (ideal) noise stability of the methods used for transmitting waves in 
terms of the intensity of the noise at the output of an ideal receiver; he has demonstrated that the expression for the 
noise intensity depends on the modulation system which is used (direct or integral modulation), 


For direct modulation systems (amplitude and phase modulation systems are the direct modulation systems 
which are used at present in radiotelemetering systems) the expression for the square of the noise intensity at the 
output of the ideal receiver is written as [1] 


On te (1) 
- ap 


If the width of the passband AF for the i-th channel is equal to the maximum frequency of variation F,,,, 
for the measured parameter and if the variation of the measured parameter is within the limits + 1, then such a 
noise intensity at the output of the ideal receiver will produce a mean-square error 


2 o 


F max 


for the measurements, 


Substituting the value Fy, = ad » we obtain 


62 = o = 3 
i ~TedAanea 2 
or | (A, oy 


’ 





on 












T/2 
I=TA,(d,t) = \ A, (h, t)dt 
—T/2 


is the specific energy of the oscillations, 


It should be noted that the resulting expression for determ ining the mean-square error of the measurements in 
RTS-FDC when direct modulation systems are used coincides exactly with the well-known expression for the mean- 
square error of measurements that are taken for transmission of instantaneous parameter values (1). 


When integral modulation systems are used (for example, FM) the expression for the noise intensity at the 
output of an ideal receiver for transmission of waves is written as (for AF =F .,,) 


din at A ee 
’ 24 (A, ay ans (A = J A(t) dt), 
OA 


The square of the effective value of the amplitude of the noise wave at the output will be 





"max (2 y 2 
\ @df a <n . 
J y 3 | OAC, ay max 
. — OA 
y under these conditions, 
Ing 


Since variation of the measured parameter occurs within the limits + 1, the expression for determining the 
square of the mean-square error for measurements in RTS-FDC when integral modulation systems are used will be 








written as 
e 
ae o? 2 max 
or [24 OP a (2) 
OA 
It is evident from this formula that the mean-square error for measurements taken in integral modulation 
) systems increases in proportion to frequency, whereas for direct modulation systems it is independent of frequency. 


Telemetering systems are often evaluated in terms of the meamsquare error that is normalized to the entire 
range of parameter values. When such a criterion is used expressions (1) and (2) are rewritten as follows: 











ia ca m3 (1") 
‘i aap 8 

sr| et | 
2 o? 2 max (2") 
i dA, 2 — 

sr | a | 


In order to obtain the values of the mean-square error for measurements that are made in a system with a 
practical receiver we make use of data on the signal-to-noise ratio (U;/U,) at the receiver output for various trans- 
mission methods cited in the book by* M, H, Nichols and L, L, Rauch [2]. 


The mean-square error of the measurements (normalized to the entire range of parameter values) shall be 
determined from ine following formula for reception by a practical receiver: 


* The data obtained by Nichols and Rauch have been verified by numerous authors and do not give rise to any doubts, 








In determining the real noise stability of any particular system it was assumed that the AM demodulator con- 
sisted of a linear detector and a filter, and that the FM or PhM demodulator consisted of a limiter, a frequency 
(phase) discriminator, and a filter, (The notation PhM denotes phase modulation.) Here it was assumed that the 
demodulator elements did not introduce additional distortion or noise, 


In this paper we have used the following notation: A(), t) is the signal which is a function of time and the 
measured parameters; AF is the width of the passband in the channel receiver; Fy is the maximum deviation of the 
sub-carrier frequency; F, is the frequency of the sub-carrier for the i-th channel; Fy,ax is the maximum frequency 
of the measured parameter; f , is the maximum deviation of the carrier frequency; M is the amplitude modulation 
coefficient of the carrier; m is the amplitude coefficient of the sub-carrier; N is the number of channels; T is the 
time interval that is associated with the maximum frequency of variation of the measured parameter by the re- 
lationship T= ¥, Fmax' Uegg 1s the effective value of the unmodulated carrier voltage; Uegs , is the effective 
value of the unmodulated sub-carrier voltage in the k-th channel; U,, is the amplitude value of the unmodulated 
carrier; ) (t) is the measured parameter that varies within the interval + 1; 6 is the meamsquare error; o is the per 
unit fluctuating noise intensity, v/VCps; uy» is the angular frequency of the carrier; wg is the maximum angular 
deviation of the carrier; 9; is the angular frequency of the sub-carrier for the i-th channel; Qy is the maximum 
angular deviation of the sub-carrier; ®, is the maximum deviation of the carrier phase, @ 4 is the maximum de- 
viation of the sub-carrier phase, The subscript "k" is used to denote quantities corresponding to one channel, 


1, The AM-AM System 





The potential (ideal) noise stability. When using double AM the signal can be written as 
A(A, t) =U,» {1 + M [(1 + mA) cos (Qit + p)]} cos (@ot + @). 


In order to find the error we determine the integral 


COS? (Wot + @) sin? (Qyt +- »p) dé 





On 


I as 24 |" Q, 2} dt = U2,M*m? 
-T/2 . 


+ T/2 
T/2 


| r2 ’ 
=F U},M?m*T . 


The value of the integral is obtained on the basis of the fact that the average value of the product of two 
functions that do not have identical frequencies is equal to the product of the average values of the factors, 


Moreover, it was assumed that 7’ S- = and 7S = . In that case it is possible to assume that [1] 


5 








cos? (Wot + @) = > , sin? (Qit + wp) = ~* , 





COS* (Wot -|- @) sin? (Qit +- ) = = 


We shall make use of these equations without reservation in our subsequent analysis. Substituting the value 
of I into (1") and assuming m =1 and U,, = V2Uegf, we obtain 


a“ re * 
eff 





The formula for 5} in the case of a multichannel system can be obtained if instead of the over-all carrier 
modulation coefficient M we substitute the value of the modulation coefficient M, which applies to one channel, 


In the general case we can write My; = 5" 
i 
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3) 








U 
Substituting the value of My into (3), making the substitution Ueg; , = rT an » and assuming approxi- 
1 


mately that M =1, we obtain 


3? 
i = <n (4) 
40 eet. 
The practical noise stability. Based on the analysis of the passage of signal plus noise through an AM-AM 
receiver, we derived the following expression for the signal-to-noise ratio at the output of the channel receiver [ 2): 


(a!) _ Vet. & 
YnJout 25 VAF 


The expression for the mean-square error of measurements made at the output of the practical receiver is 
written as 





5 (2) i (5) 
8 |- ig 
p U, 20 tf k 


Using I], to denote the error ratio between the ideal(Kotel‘nikov) and practical receivers and taking into 
account the fact that T = */, F_.,,, we obtain 





yr — °> — _4F 
2 F 
1 max 


From the derived formula it is evident that for AF =F max the practical noise stability of an AM- AM system 
reaches the ideal noise stability. 


2, The FM-AM System 





The ideal noise stability. When FM-AM modulation is used the signal can be written as 


A (hd, t) = Um {1 -+- M cos [(Qyt + QgA) 4- Py)}} COS (Wot + Go), 
where A = {A (f) dt. 
The per unit energy of the wave A’, (\, t) will be 


ee. +T, 
ae < (24 G OP gy . Um 6) 
oa |] ~T \ LOA od. 


Substituting the derived value into the formula (2*) and making the substitution U,, = ¥ 2Ue¢¢, we obtain 


§? - 1 3 


g OF ee eee Qinax (1) 
i 12 UZ M*QeT 


In order to obtain 57 for a multichannel system it is necessary to make the same transformations as those 
made in section 1. 


Then 


& aN o? oy 
Nias “42 U2 QT max’ (8) 
d 
e 
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The practical noise stability, The expression for the meamsquare error of the measurements at the output of 
the practical receiver is determined by the expression [ 2) 


2  1,;Un\? 1 AQ%*#AF 
%=3(00) =o os m 
: eff.k 4 


Using I, to denote the ratio between the errors for the ideal (Kotel'nikov) and practical receivers and taking 


into account the fact that T =} Finax? We obtain 


. AF? 


qe = 


alate 





3 
* max 


| 
- t 


Thus the practical noise stability of an FM-AM system is equal to the ideal noise stability for AF = F_. .. 
3. The PhM-AM System 





The ideal noise stability. In order to find the mean-square error we write the signal in the form 


A(A, t) =U, {1 + M cos [(Q)t + 8 ,A (t)]} cos (Wot + ¢), 
where 





Determining the integral I and substituting its value into (1') for U,, = ¥2Ueff , we obtain 


. 32 32 
Dt ee (10) 
i 8] 2742 92 ‘ 
ge oT 
Performing the same transformations as those performed in section 1, we obtain 
2 
eo = ———. (11) 
. a. ae 
eff.k 4 


The practical noise stability, The expression for the mean-square error of the practical receiver is written as 


& — otAF 


Ts 2 2 ° (12) 
20 © a 


The error ratio for T = */, F,,4, will be 


m — > _ AF 
of F max 





Thus, the practical noise stability of a PhM- AM system for AF =F ax Wiil be equal to the idea] noise stability. 
4, The FM-EM System 
The ideal noise stability. We shall represent the signal in the form 





A (A, t) = Um 608 [(Wot + Q4t +- QA) + Pol. 
The per unit energy of the wave A‘, (\, t) will be equal to 
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Substituting this quantity into the expression for the error (2") while taking into account the fact that Uppy = 
= ¥ 2Uefe k » We have 





#Q? 
6 = Fr (a3) 
1 94 
24 UN at 


The practical noise stability. The signal-to-noise ratio at the output of the channel receiver of such a system 
can be written as 


U; j a0 Qq 
t= (4) 
wT. eff.kA2 Vo5VAF | 


The mean-square error for the specified method of reception will be equal to 


A (aps AFe 


12 u2 QQ? 
eff.k @ 


6, = (15) 


The error ratio for ideal and practical reception is 


5, A F3 
I? = ei. oe AF : 
oF 3 
max 


Thus, this type of radio-receiver also realizes the ideal noise stability in practice. 
5. The PhM-EM System 





The ideal noise stability. The signal is written as 
1(A, t) = Um, C08 [(@p 4- Q))t 4- Q qa’ (t)). 


The mean-square error is determined (just as in section 4) by the expression 


oF = i a ° (16) 
ert d? 
The practical noise stability. The expression for the signal-to-noise ratio at the output of the practical 
receiver is written as 
(-3 ) Ue peen 
\ 7] dy - 
Un: out eff.k V¥23VAF° 


The mean-square error of the measurements will be determined (for AF =F,,4,) by the formula 


3F my 
5;, - =. (17) 
4U? §\- 
eff.k d 


The ratio between the errors corresponding to practical and ideal reception is 


Thus, a practical receiver in a PhM-EM system also realizes ideal (for AF = F,,,,) noise stability of the tans- 
mission method, 


6. The AM-FM System 
The ideal noise stability. For this system the signal can be written as 
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A(A, t) = U,, cos J enol +. w4)\ {1 + mA (t)} cos (Qit + wp) dt + 9} ‘ 


Usually the maximum frequency of variation of the measured parameter is much less than the frequency of 
the sub- carrier (i.e., a. « Qj). For such an assumption the expression for the signal in an AM-FM system can 
be rewritten in the following form after taking the integral: 


A(A, t) = U,, cos { 0 -t. = [1 + mA (t)] sin (Q,¢ + p) + 9} : 


The mean-square error of the measurements for such a signal can easily be found according to formula (1), 
For m =1 and U,, = V 2U ese. 


83 3 Q? 
Bee 4Ufpo*,T ° (18) 
The expression for 5] in a multichannel system can be derived from expression (18) if instead of the quantity 

wg we substitute the value of deviation corresponding to a single channel (i.e., wg ,) and replace Uggs by Uege x: 


3? 0? 


2 2 " 
oy? 


2 
hater 
éff.k 


(19) 
The practical noise stability. The mean-square error of the measurements performed for reception of AM-FM 


signals with a practical receiver was derived on the basis of the data cited by Nichols and Rauch [2] and is deter- 
mined from the expression 


32 o* OF AF 
Pe 372 Way 
Ure Od 


The ratio of the errors for reception by practical and ideal receivers (for T =*/, F,,,,) will be 


ae 


ae ae 


Il 
ge OF 


cr 


max 


Thus for AF =F pax the practical noise stability of an AM-FM system will equal the ideal noise stability. 
7. The FM-FM System 
The ideal noise stability. The signal for such a system can be written as 





A (A, t) = U,, cos {oot + «4 \ cos (2 + QA +») dt + 9} - 


‘al \? (t) dt. 
It is very difficult to find the expression for the mean-square error in general form for reception of such a 


signal by an ideal receiver. Therefore, we shall first simplify the expression for the signal. 
Assuming that Q; > 0,,,, and taking the integral, we obtain the following expression for the signal: 


A (i, t) = Um 008 {ot + isin (Q;t + Q,A +) +q}. 


The mean-square error of the measurements for reception of such a signal can be determined from formula (2). 


1400 





ty 


)) 





6a 


The per unit energy of the wave A‘, (), t) will be 


+T/2 





1 2 ee 
T ; OA 49? 
-—T/2 


Substituting the derived value into the formula for the error of a multichannel system, we have 





20202 
82 mf G 25 2" ax (21) 
i 12 yy? ot Q?T” 
eff.k @ 4d 


The practical noise stability, The expression for the square of the mean-square error for reception by a 
practical receiver is written as 


o 3604 QF (2x)? AF* 
0% =F 2 toe 
Yip n %a% 


The error ratio determined from formulas (21) and (22) will be 


(22) 





. é? 3 
i max 


From the derived relationship it is evident that for AF =F,,,, the noise stability of the practical receiver 
reaches the ideal nose stability of the transmission method, 


8. The PhM-FM System 





The ideal noise stability. The signal is represented in the following general form for such a system: 


A (A, t) = U», cos {toot > w4\ cos [(Q;t + 04d (t) + p] dt + a} - 


Since usually the frequency 9; of the sub-carrier is much greater than the maximum frequency Q,,,, con- 
tained in the message } (t), it follows that the analytical expression for the signal can be simplified if we write 


A(h, t) =U, cos {oot 4 a sin [Q;t +0, %(t) + p] + 9} (23) 
a 
The mean-square error of the measurements for such a signal can be found from formula (1), The expression 
for the error in a multichannel system will be written as 
65 = (24) 
© Wee eegOa? 


The practical noise stability. On the basis of the data in [2] the expression for the square of the mean-square 
error can be derived in the following form: 





of Q? AF 
2 i 
5p os —a «6Laa . (25) 
eff.k 


The ratio between the squares of the errors will be 


of AF 
f——j pe = —_— 
M=3 
i max 
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For AF =F ax 4 practical receiver in a PhM-FM system realizes the ideal noise stability of the transmission 
method, 


9. The AM-PhM System 





The ideal noise stability. The signal is written as 


A(A, t) = Up, cos {@ot + Da[(l + mA) cos (Q;t + p)}}, 
where O4 = w4/Q,. 


We shall determine the integral 


+T /2 


I = 7 [Ae at = Umea, 


The mean-square error for m=1 and U,,, = V 2Ue¢, is 


pe (28) 
i 8! on?o27T 
te 


In order to obtain 6; for a multichannel system we must replace }q with the value oq , (i.e., with the value 
of the phase deviation for one channel), In the general case $4 ,=q/ f(N). Assuming Uses 1, = Ucss/fs(N), 
we obtain 


2 o 
Oe ee 5 
1 4u? oT (29) 
eff.k 


The practical noise stability. The square of the mean-square error of measurements taken for reception with 
a practical receiver in an AM-PhM system is determined from the expression [ 2] 





U._\? 2 

ret 7 ay —- oAF 

%=% (7 }. 7 eT ee (30) 
eff.k 4 
The error ratio will be 

& AF 

Tr = oak = e 

of F max 


It is evident from this formula that for AF =F ,,,, the practical noise stability of an AM-PhM system is equal 
to the ideal noise stability. 


10. The FM-PhM System 





The ideal noise stability. The signal is written as 


A (A, t) = Up, cos {wot + Dy cos (Qt + Q yA + p) + 9}. 


Based on formula (2°) the expression for the error in a multichannel system is written as 


: (30°) 
i {2y? @98T 


The practical noise stability. The expression for the square of the mean-square error for practical reception 
of FM-PhM signals is written as [2] 
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ue 











55, so (31) 
6 772 2 m2 
Cte.Keo Ma 
The error ratio will be 
Bie oF . AP 
7 5 ae ae 
. P’max 
For AF =F ya, the practical noise stability of the FM-PhM system reaches the ideal noise stability, 
11. A PhM-PhM System 
The ideal noise stability. The signal is represented in the form 
A (A, t) = Um cos {Wot + D ,cos [Qit + 8,4 (t)]}. 
Based on the formula (1') we find the mean-square error in a multichannel system: 
2 
i U2 @*, 92 T 
: eff 4 a 


The practical noise stability. The expression for the mean-square error of the measurements of the output 
from the practical receiver in the system is written as [2] 


2 oAF 


4, = = = 
F 9172 2 
Vite d® d 


(33) 


The ratio of the errors for practical and ideal reception will be equal to 


Thus the practical noise stability of a PhM-PhM system for AF =F,,,, 1s equal to the ideal nolse stability. 


SUMMARY 


The paper compares various types of multichannel radiotelemetering systems with frequency division of the 
channels in terms of ideal and practical noise stability, The results of the comparison demonstrate that the practical 
noise stability of the radiotelemetering systems under study may reach the ideal noise stability when weak fluctuating 
noise is present. One of the necessary conditions for the practical noise stability to approach the ideal noise stability 
is the use of an ideal low-pass filter with a cut-off frequency equal to the maximum frequency of the measured 
parameter (AF =F,,,,) at the output of each channel; i,e,, optimal filtering must be assured, Since practical filters 
never have an infinitely steep cut-off, it follows that the practical noise stability of a system will always be some- 
what worse than the ideal noise stability. 


The author expresses his deep appreciation to A, P. Manovtsev for a number of valuable remarks which were 
made during the preparation of this paper. 
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A method is proposed for the investigation of automatic-control systems, based on the 
transformation of the coefficient matrix of the differential equations describing the system in 
question, The method can be used to study the stability and quality of a system by simple 
involution of the transformation matrix, 


We consider the system of linear differential equations 


dx 
—— = anti + Gite +... - + 4yyTp, 
at 
dx 
——* = agi%1 + Ggere + ..- + Agntn, 
dt 
oo . . . * . . . . . . .* (1) 
dx ’ 
= = QyyT1 1 Angt2 + --- T @ea%n, 


where t is the independent variable, x,, Xz, ..., Xp are the functions of the independent variable t , and aj 
(i, j=1, 2,..., m) are real numbers, 


This system can be written in the form of a single matrix differential equation: 


dx 


at = Az, (2) 


where A is the matrix of the coefficients of the system (1) and x is the column matrix with elements x, X, . . . »Xp- 

The matrix A, corresponding to the system (1), can be tansformed by a conformal mapping of the form 
where R is the radius of a circle containing the eigenvalues ); of the matrix A, m=R-+ a is the coordinate of the 
center of the circle, and a is the assigned degree of stability. 


The transformation (3) maps a circle of radius R in the left half-plane of the complex variable ) into a unit 
circle in the p-plane, 


We note that various conformal transformations of the plane of the complex variable ) are possible, On the 
- unit circle we can map the whole left half-plane of ) [1], or a part of this half-plane bounded by rays, curves of 
the second-order, or other curves, These mappings lead to various transformations of the matrix A, and naturally 
to various magnitudes and qualities of the results of calculations, In what follows, we will restrict our considerations 
to the transformation (3), as the best for application in practice, When ) from the formula (3) is substituted in the 
characteristic equation | A— AE | =0 of the original system of differential equations, the result is 
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A 
A=Rt+E and A=A+aE, 


When the solution of the system (1) is asymptotically stable, then * absolute values of all the eigenvalues 
of the matrix A, will be smaller than one (| p;| < 1) and we will have aN + 0 for Noo, Thus the successive 
powers of A, can be used to determine the stability of the system, and there is asimple region of stability. Here 
the stability for the Nth power of the matrix x As is calculated by means of any of the norms of the matrix (\|A"\\|<‘1) 
and the instability from the criterion |S pAN | > n, where n is the order of the matrix, 


The tangent to the circle at the origin in the )-plane can be used to calculate the oscillation of the system, 
where the oscillation index is a function of the ratio R/a, i.e., each pair of values of R and a corresponds to a 
definite oscillation of the system. 


It should be remarked that the calculation of the matrix A, does not require the inversion of the original 
matrix, The matrix A, is obtained from A by elementary arithmetic operations, 


The obtaining of better qualitative results, and especially the obtaining of a complete picture of the quality 
(curves for transient processes) of the given system, obviously depends on how accurately our transformation matrix 
reflects the structure and properties of the integral matrix of the system of differential equations being studied, 
The effectiveness of various transformations of the original matrix, among which is the transformation of the type 
(3), can therefore always be evaluated by comparing the transformed matrix with the canonical structure of the 
integral matrix of the system being investigated, The integral matrix of the Eq, (2) is the matrix of a certain 
fundamental system of solutions of the original Eq. (1) (2). 


We denote it by X(t): 


wir (t) win(t)..+ 2p (t) 
wor (t) woo (t)..- Zon, (t) 


4.8.4: 024076 # © &£.5 3 9.9. € ae 


Tn, (t) Inq (t)--+ Tan (t) 


X () = 


In the case of a system of the type (1), the matrix A commutes with its integral [3], and a family of integral 
matrices, satisfying Eq. (2), can be given in the form 


X (t) = e4' Xo, (5) 


where eAt is the integral matrix, and X, is an arbitrary, constant, nonsingular matrix, 
If X, is the matrix of initial values, then formula (5) gives the general solution of the system (1), 
The logarithm of the matrix A, can be expanded in a series: 





= 
In Ai = > f 2. (Ai — E)*, (6) 


k=1 
since §=[A,—1/<1 (i =1,2,...,8) 
For simplicity, we set a = 0, and then 


e—-L< +...) 
AN = eNinAs eg \RO 2 RTs, (7) 


If we consider only the first approximation, then we find that the matrices eAt and AN have the same structure, 
The difference lies only in the fact that the role of the time is played by the discrete parameter N/R, Thus by 
successively raising the matrix A, to higher and higher powers and operating with the resulting matrices on the vector 
of the initial values, we can obtain numerical solutions of the original differential equations (curves for transient 
processes), Here the value of N/R corresponds to the time interval. 
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The relative error of the calculations increases with increasing powers of the matrix (or with increasing time 
t=N/R), The choice of R must ensure the prescribed accuracy of calculation, within the limits of the necessary 
time of integration, The calculation is easily carried out using the norm of the matrix. Using the first terms of 
the expansion (7), we obtain 


uy 
oR” , (8) 
net CO) 


N 


where § is the given relative error of the calculations, 


Thus the maximum attainable time for which the given accuracy of integration is maintained is determined 
by the relation 





OR 
‘max °' [aay (9) 
where 6; = In i 
— 


Since the norm often leads to results that are not accurate enough, it is useful to transform the calculation to 
one directly in terms of the eigenvalues of the matrix: 


a1 x 4%? 
inpp=in(i+st)— 4-4... (10) 
\A max? 
a 1 (11) 
e Ses ° 


2R 
t max 5 7x’ (12) 
max 


where R,,,;,, is the minimum radius for which the circle still contains the eigenvalues ) ;. 


The formulas (9) and (12) can also be used to obtain the value of the maximum allowable power of the matrix 
for which the piven accuracy of calculation is still maintained, 


Both for the stability and for the analysis of the quality of the system, the powers of the matrix must 
follow the law N = 2% (k=1, 2, 3,...). The selection of the value of R by using an electronic computer presents 
no difficulties, Large values of the radius increase computing time, but the number of steps that must be carried 
out by the machine increases only at the rate log,R. For maintaining the prescribed accuracy for large values of 
the time, it is convenient to double the integration grid, 


Practical testing shows that the method is very effective in investigating automatic systems described by 
differential equations of high order, 


In conclusion, the author wishes to express his sincere thanks to V. I. Zubov and N. S. Levchen for posing the 
problem, and to F, R. Gantmakher for his valuable comments, which were taken into account in preparing this article, 
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ELECTROMAGNETIC CONTROL ELEMENTS 


I, M. Krassov and V, N, Nikol'skii 

(Moscow, Kuibyshev) 

Translated from Aytomatika t Telemekhanika, Vol, 22, No, 11, 
pp. 1546-1549, November, 1961 

Original article submitted November 9, 1960 


We discuss the results obtained in the development of new electromagnetic control 
elements with progressive armature motion, 


In modern automatic regulatory and control systems electromagnetic devices are widely used as control 
elements, 


As a result of the increasingly severe requirements which automatic-control systems must satisfy,special 
importance has been attached to the static and dynamic properties of electromagnetic devices, The design and 
manufacture of improved types of electromagnetic devices has become a necessity, 


G52 








6 


Fig. 1. Construction of electromagnetic control element, 

1 and 2) Coils of the first and second electromagnetic systems; 
3) armature; 4) magnetic conductor (core); 5 and 6) balanc- 
ing springs, 


In this paper we present a brief description of the results of a study aimed toward the development and in- 
vestigation of new electromagnetic control elements, The problem was to try and develop electromagnetic control 


elements with improved static and dynamic properties which were lightweight, of small dimensions and could 
operate under vibrational conditions, 


The control elements which were developed consist of electromagnets (Fig. 1) with a two-sided linearly-dis- 
placed armature, Two electromagnetic systems act upon this rotor in opposite directions, The armature is held in 
a central (neutral) position by means of two springs, For all other armature positions the total force, which must 
equal the torque of the electromagnet,represents the sum of the spring forces and the external load, 


1407 








The photograph (Fig. 2) shows the external view of the electromagnets. The magnetic circuit of the electro- 
magnets is made of Armco steel. In order to decrease the eddy currents in the core there are slits in the cores, the 
covers and the armature, The balancing springs are made of phosphor bronze and are isolated from the covers by 


means of textolite rings. 





Fig. 2, 


The usual technical conditions which the electromagnets must satisfy consist of a series of to a certain extent 
contradictory requirements, such as, for example, the requirement of a light weight and large torque, the require- 
ment of low power usage and small time constant, etc, Therefore, the choice of the type of electromagnet and its 
basic dimensions is determined by the principal requirements and its details are corrected according to the remaining 
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Fig. 3, Static characteristics of one of 
the electromagnetic control elements. 
Coil currents I, and I, shown as per- 
centages of the maximum current. 
Armature displacement x as a percent 
of the maximum displacement. 
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requirements, In our case we took as the basic requirements: linear 
torque characteristics and a high self-oscillation frequency; these char- 
acterize the dynamic properties. 


In order to satisfy the first requirement we chose an electromagnet 
with a cylindrical armature and conical ends and a stationary core with 
conical poles. Investigations have shown that this shape of pole enables 
us to obtain linear torque characteristics for relatively large torques. 


In order to satisfy the second requirement, i.e,, in order to obtain 
a sufficiently high self-oscillation frequency we must decrease the weight 
of the moving system of the electromagnet. Some of the electromagnets 
available at the present time, in particular [1] and [2] have movable 
systems consisting of two armatures with a common axis and stationary 
systems in the form of two nonmagnetically-associated electromagnetic 
systems. In contrast to this type of construction,in the magnets which 
we have developed both electromagnetic systems are connected by means 
of a common pole with an inactive (parasitic) air gap and have a common 
armature, This permits us to decrease significantly the weight and di- 
mensions of the movable part and also the weight of the entire electro- 
magnet, 


For the rotor displacement to be a linear function of the control 
signal for a chosen form of the electromagnet armature and poles,the 
magnetic reluctance of the parasitic air gap was reduced to a minimum 
by proper choice of air gap dimensions, For this condition and for low 
saturation of the steel core the forces acting upon the armature are given 
by the formulas 


dG 
Fy =c (Iw, (1) 


Fy =e (Iw) 52, (2) 








ing 


et 
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ets 
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where F, and F, are the forces arising in the first and second magnetic systems, I, and I, are the currents in the coils 
associated with these systems, w is the number of turns in each coil, G, and G, are the magnetic conductivities of 
the operating air gaps of the system, x is the armature displacement from the mean position and c is a constant 
coefficient. 


For oscillatory operation the electromagnetic control element operates together with a vacuum tube amplifier. 
Here constant voltage pulses of the same or varying duration are applied alternately to the electromagnet coils. As 
a result the armature will oscillate about a position determined by the pulse duration and dynamic properties of the 
electromagnet. Since the pulse duration is a function of the input control signal, by changing the input control 
signal we can change the values of I, and 1, proportionately and therefore the displacement of the oscillating 
armature with respect to the neutral position. 


When the control signal is changed the currents in both coils change, the current in one coil increasing and 
that in the other coil decreasing so that their sum remains unchanged. Therefore, we can write; 
h=1,,", (3) 
I,=1,,(—n), (4) 


where In is the maximum value of the current in each coil and a is a coefficient: 0 <n <1. 


The resultant force acting upon the armature is the difference of the forces obtained from the two systems 


dG dG 
F = F,—F,=c[(1—n)/,,v]*  —c [nl,,w}* (5) 
With the additional condition 

dG, _ dGs 

dz ~ dz (6) 
the expression for the resulting force simplifies to 
dG 

F =e (/,,w) (1 — 2a). (7) 


Formula (7) shows that the force is a linear function of n and therefore a linear function of the control signal. 


Condition (6) is fulfilled for the allowed armature displacement in the electromagnets which we have de- 
veloped by the choice of the form of the armature poles and core and by the large size of the initial air gap. 


The derivatives dG,/dx and dG,/dx are determined from the expression for the conductivity of the operating 
84p: 


G 





___ 1,256. 2ar m ; 
1.2 (8 + 2) cosa | sina — (+ 2) sina], 


(8) 


where §, is the initial operating gap corresponding to the neutral position of the armature and measured from the 
armature axis, a is the angle at the base of the armature cone, m is the height of the armature cone andr is the 
armature radius at the inner cut of the core pole. 


If we expand the expression for the resultant force F, obtained by taking formula (8) into account, in a Taylor 
series and then discard all terms involving powers of the variable higher than the first, we can write 


F = k,Al + k,Az. (9) 


Here k; and k, are constant coefficients determined by the construction of the electromagnet, AI is the change 
in coil current and Ax is the change in the armature position. 


1409 








The equation for the spring force N can be written 

























N=c,2, (10) 
where the c,, are the coefficients of stiffness of the springs. Equations (9) and (10) permit us to cotain a linear 
relationship between the armature displacement and the current changes in the coils. 


In Fig. 3 we show the static characteristics of one of the electro- 
magnets; a plot of the armature displacement versus the coil currents. 

















I—— <— The curve is based upon the experimental data and confirms the linearity 
5 \ of the function. 


Figure 4 is an oscillogram of the armature displacement from the 
neutral position to the boundary position during a discontinuous change 
in the voltage applied to the coils. This oscillogram characterizes the 
dynamic properties of the electromagnet. 


Fig. 4. Oscillogram of the transient 
process. 1) Armature displacement; 
2) changes in current; 3) time 
(per 1/500 sec). 

As an example we give the basic technical data for one of the 


electromagnets that we have developed: 





total weight of assembled electromagnet 260 g, 
weight of armature and rod 8 g, 
armature displacement within the limits of the 

linear characteristics +1.3 mm, 
torque for a 1.3 mm armature displacement 1000 g, 
power used by electromagnet 5.7 wW, 
self-oscillation frequency 120 cycles, 
time required for armature displacement 

from neutral to boundary position 0.032 sec. 
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cover English translations appears at the back of this issue. 














EVENTS 


SCIENTIFIC SEMINAR ON PNEUMO-HYDRAULIC AUTOMATION 


The Fourth All-Union Conference Seminar on Pneumo-Hydraulic Automation was held at the Institute of 
Automation and Remote Control. 


321 scientists and engineers and 97 organizations working in the field of pneumo-hydraulic automation in the 
various cities of the Soviet Union attended the sessions. In all,representatives from 19 cities attended the conference, 


In his opening address, Professor M, A. Aizerman, the seminar chairman, briefly listed the successful new 
developments in pneumo- and hydro- automation in the Soviet Union during the last few years, It was shown during 
the First International Congress, IFAK,on Automatic Control held in Moscow last summer that in the field of pneumo- 
and hydro- automation we are not only keeping abreast of the foreign scientists but are ahead of them as far as some 
questions are concerned. 


Professor Aizerman furtherremarked that two volumes of papers were published during the operating period of 
the seminar. The seminar provided a great impetus to the development of pneumo-automation and hydro- automation, 


After the opening address A. A. Tal* (Moscow) and M, I, Zhutovskii (Moscow) told of a new system of elements 
used in industrial pneumo- automation, Developed by the Institute of Automation and Remote Control of the Academy 
of Sciences of the USSR and presently being used by the “Tizpribor* factory, the universal system consists of ten 
standard elements; the system permits us to set up the various automatic-control circuits involved in the different 
production processes, 


All the elements have small dimensions (30 x 50 x 50 mm on the average). The location of the element 
outputs is unified so as to permit the use of standardized, uniform components. The units of the assembly, which 
consists of a group of functionally-related elements, are mounted on special flat plates with internal communication 
channels. 


The paper presented by V. P.Temnog (Moscow) described the work that he and V. A. Khokhlov (Moscow) had 
carried out jointly in the design of a hydraulic regulating column (HRK-1) for general industrial use. 


The paper outlined the principles underlying the design of complex hydraulic regulating systems with electric, 
pneumatic or hydraulic input signals. Due to the special feed methods used, the column, which applies a powerful 
hydraulic drive, requires only 100 watts of energy. 


N. P. Zhivov and Naizer (Moscow) described a pneumatic information machine which was composed of 
pneumatic logic elements, 


In their paper, “Membrane Type Computing System, Operating at Low Pressures (0-100 mm of water) G. T. 
Berezovets and V. N. Dmitriev (Moscow) showed that by using a combination of operational amplifiers, constant 
and regulating valves and low pressure pneumatic chambers they could construct control and regulatory devices 
and also computers which can perform algebraic calculations and can differentiate and integrate. 


Sh, I, Israelov presented a paper written jointly by A. A. Abdullaev, K. A. Oganova and Sh. I, Israelov 
(Sumgait) entitled, "Pneumatic Calculator and Drive for the Automatic Control of the Multiplicity of Circulation 
in Processes Involving the Circulation of Finely- Dispersed Catalysts.* 


Professor Israelov described the elements that had been developed for the automatic regulation of a complex 
parameter of the multiplicity of circulation of the catalyst during the dehydration of butane in a “boiling” layer of the 
circulating finely-dispersed catalyst. 


In his address Sh. 1. Istaelov described the computer which had been developed by the authors. It consists of 
two multiplier-dividers and was developed at the IAT of the Academy of Sciences of the USSR, He also described 
4 specially-constructed valve with an automatic drive. 


An automatic system that controls the multiplicity of circulation has been developed using the above-described 
elements, As a part of this system there has been developed a system that automatically controls the level of the 
Catalyst in the reactor. 











V. V. Aronovich (Moscow) described the application of pneumatic logical elements and extremum regulators 
in automation schemes for use in the chemical industry. 


V. P. Zenchenko, in his paper “Receiving and Actuating Devices in Pneumatic Control Systems for Automatic 
Machinery," examined the general problem of the construction of a circuit for the control of automatic machines 
and indicated the method of solution of the problem. In his paper he studied the universal pneumatic push-button 
devices, If we include several of these push-buttons we can perform several logical operations. He described a 
pneumatic drive, a piston whose rod could occupy various positions depending upon the state of the pneumatic circuit, 


Ts. P. Krivol (Moscow) presented a paper, “Pneumo- Automatic Controls for Home Martensite Ovens," 
Together with TsPKB and the factories, where a trial automation will be carried out, the TsLA has developed auto- 
matic controls which will distribute the air blasts (air stream) along the pipes and forms used in domestic ovens and 
which will also provide automatic regulation of the heating conditions in martensite ovens, The circuit (system) 
for automatic control of the air stream in the air shafts of the home oven was developed at the metallurgical factory, 


B, S. Tarkhovskii and L, O, Khvilevitski (Moscow) in their article "A Pneumatic Correlator" describe a 
machine designed to process low-frequency information. This machine is to be used not as a laboratory instrument 
but rather as a control element which would be part of an industrial control system. The application of pneumatic 
computer elements permitted Tarkhovskii and Khvilevitski to design a simple and successful correlator which 
satisfies contemporary requirements. 


In their communication, N. G. Luk"yanova and V. M. Eigenbrota (Moscow) examine some questions arising in 
the design and construction of pneumatic multiple- point controls and regulators, The article indicates that the 
method of dynamic compensation may be used for giving remote directions and for remote control of a regulator. 
It may also be applied in digital computers. 


Long pneumatic lines and methods of design of these lines were described by L. A. Zalmanzon (Moscow). 


E. V. Gerts and G. V. Kreinin (Moscow) in their paper, “Design of Pneumatic Drives" outlined a method of 
design which facilitated the choice of dimensions for the elements of the machine and the calculation of the 
operating cycle of the pneumatic system. 


The talk of Yu. N. Gerulaitis and I. V. Serdyukova (Moscow) was devoted to a membraneless transducer de- 
veloped at the NIOPIK for transforming small forces into compressed air pressures, 


A number of examples were given of the use of the membraneless transducer. The simplest of these was the 
application of the transducer to obtaining a pneumatic signal at the output of a scale (in the construction of weight- 
measuring machines), Another example of the application of a membraneless transducer was a measuring technique 
where an electronic- pneumatic circuit was used as a transducer, other examples given were magnetic gas analyzers 


for hydrogen, etc, 


Ya, M, Mar*yanovsky (Moscow) spoke on electropneumatic transducers. In order to build a transducer which 
did not require an amplified input signal one must increase the sensitivity of the compensating devices. One solution 
of this problem consists of using the dynamic action of a stream of air from a nozzle as part of a feedback loop, 


The paper of E. A. Andreev and L. A. Tenenbaum (Moscow) was devoted to the pneumatic transducer, The 
transducer used the force of an air stream flowing from a nozzle onto a slide or shutter. The transducer consists of 
two nozzles (low and high pressure), a self-centering shutter (or slide) which is suspended in the air stream and a 


high-pressure pneumatic valve. 


E, M, Freedman (Kiev) described in his paper experimental investigations of high-temperature pneumatic 
pistons for hydraulic drives and pneumatic followers. The results of his investigations showed that the output and 
load characteristics of the piston in the hydraulic drive studied improved significantly with temperature increases, 


Yu. P. Zolkina (Moscow) described a switch which made use of a pneumatic relay developed at the Institute 
of Automation and Remote Control of the Academy of Sciences of the USSR. The switch was to be used in an 
automatic test system which forwarded samples obtained at several test points to the input of a gas analyzer. 


M. L, Podgoetski (Moscow) presented the results of a study of the circuit of a compensating pneumatic trans- 
ducer with an isodrome amplifier. 
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In the section on hydro- automation the first paper read was entitled, “On the Stability of a Regulatory System 
Controlled by a Hydraulic Servomotors with Stops". It was presented by Professor A, M, Letova. 


G. N. Knyazev (Moscow) delivered a paper on "A Study of the Stability of a Hydraulic Drive with Strong 
Feedback.” 


V. N. Baranov (Moscow) and Yu. E, Zakharov (Kaluga) examined the forced oscillations of a hydraulic servo- 
motor with a controlled cut-off valve and strong feedback. On the basis of theoretical solutions, an engineering 
method of designing servomechanisms used in connection with oscillatory operation was developed. The calculations 
used in the design and construction of a new type of oscillator — an electro-hydraulic oscillator used as a follower— 
were presented, 


V. A. Leshchenko (Moscow) presented some of the results of a study of a hydraulic supply (hydraulic follower) 
which took nonlinearity into account. 


V. A. Khokhlova (Moscow) presented a paper on the experimental study of the volume durability of mineral 
oil used in a throttle (valve) controlled hydraulic testing mechanism. An explosion in the liquid destroys its constant 
flow, leads to the appearance of bubbles and in addition significantly increases its compressibility, In connection 
with this the durability (strength) of a liquid has an especially important meaning in the building (design) of rapid- 
acting high-quality followers. 


Professor E. M, Khaimovich (Kiev) described the work in hydro- automation of the department of jigs and 
fixtures (KIP), 


A. P, Kirpichey (Leningrad) described the results of the study of the instability of cylindrical valves in hydraulic 
systems and methods of reducing the instability by means of hydraulic centering. He also described the study that 
was made of the coefficient of leakage of oil through the working openings in the valves and negative feedback 
circuits of hydraulic automatic-control systems. _ 


The results of the study made of the hydraulic feed in transportation devices in automatic lines were presented 
by L. S. Bron (Moscow). In order to be able to study the types of electro-hydraulic regulators a study was made of 
the dynamic characteristics of the electro-hydraulic regulator of the “Tizpribor* factory which had a hydraulic 
amplifier with a stream flow. 


In his paper M. A, Yastrebenetzki (Moscow) told of an investigation which proved that in comparison with 
the hydraulic regulators previously designed at the "Teploavtomat" factory electro-hydraulic regulators have a series 
of advantages which can be exploited and also a wider region of application. The range of parameter changes in 
electro-hydraulic regulators is significantly wider than in hydraulic regulators, 


Two addresses were presented at the concluding plenary session. V. M. Gorokhov (Kharkov) delivered an 
address on a new electronic-hydraulic automatic-control system. One of the methods of perfecting automatic re- 
gulators was a judicious combination of electronics and hydraulics or pneumatics, 


The basis of the combined automatic-control system which was developed at the *Teploavtomat” factory 
was the use of electronic amplifiers and the development of the law of regulation and control for the use of hydraulic 
drives and linear electric circuits. 


Professor M, A. Aizerman described a new region of pneumo-automatics — stream techniques — and he 
described new principles underlying the design that could be applied in the development of pneumatic stream de- 
vices, which present new possibilities in pneumo-automatic devices. 


In conclusion the members of the seminar noted the significance of the papers delivered every year at the 
conference in facilitating the scientific contacts between the various scientific organizations working in the fields 
of pneumo- and hydro- automatic devices. 


A. |. Semikova 
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